//--------------------------------------------------------------------------
// File and Version Information:
// 	$Id: GqaSemileptonic.cc,v 1.47 2001/02/03 05:18:11 abi Exp $
//
// Description:
//      Module to work in the BaBar software framework for
//      analysing Knunu truth information to study the performance 
//      of the generators.
//
// Environment:
//	Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//      Anders Ryd                    Original Author
//
// Copyright Information:
//	Copyright (C)
//
//------------------------------------------------------------------------
#include "BaBar/BaBar.hh"

//-------------
// C Headers --
//-------------
#include <assert.h>

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

//-----------------------
// This Class's Header --
//-----------------------

#include "GeneratorsQA/GqaSemileptonic.hh"

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "AbsEnv/AbsEnv.hh"
#include "CLHEP/Alist/ConstAList.h"
#include "CLHEP/Alist/ConstAIterator.h"
#include "CLHEP/Alist/AList.h"
#include "CLHEP/Alist/AIterator.h"
// #include "CLHEP/Random/RandFlat.h"
#include "HepTuple/TupleManager.h"
#include "HepTuple/Tuple.h"
#include "HepTuple/Histogram.h"
#include "HepTuple/HTValOrderedVector.h"
#include "HepTuple/HTValVector.h"
#include "AbsEvent/AbsEvent.hh"
#include "AbsEvent/getTmpAList.hh"
#include "ProxyDict/IfdStrKey.hh"
#include "ProxyDict/IfdKey.hh"
#include "ProxyDict/Ifd.hh"
#include "AbsEnv/AbsEnv.hh"
#include "GenEnv/GenEnv.hh"
#include "PDT/Pdt.hh"
#include "PDT/PdtEntry.hh"
#include "PDT/PdtLund.hh"
#include "StdHepData/StdHepTrk.hh"

#include "PepEnv/PepEnv.hh"
#include "PepEnvData/PepBeams.hh"
#include "PepData/PepCollision.hh"

#include "Beta/EventInfo.hh"
#include "Beta/BtaCandidate.hh"
#include "Beta/BtaMcTruth.hh"

#include "ErrLogger/ErrLog.hh"
#include "BaBar/Constants.hh"
#include "EffTable/EffTable.hh"
#include "BetaPid/PidWeight.hh"

#include <string>
#include <map>
using std::cout;
using std::endl;
using std::ostream;
using std::map;
// using std::istrstream;
// using std::setprecision;
// using std::setw;
// using std::strstream;

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

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

// in general, a module constructor should not do much.  The beginJob or
// beginRun members are better places to put initialization
GqaSemileptonic::GqaSemileptonic(const char* const theName, 
			     const char* const theDescription ) : 
  AppModule( theName, theDescription ),
  _eventInfoList("eventInfoList",this,"Default"),
  _btaTruthList("truthCandidates",this,"MCTruth"),
  _applyBoost("applyBoost", this, true),
  _pepCollisionKey("pepCollisionKey", this, "Default"),
  _makeNtuple("makeNtuple", this, false),
  _numPrintEvents("numPrintEvents", this, 10),
  _booster(),
  _evtNum(0)
{
  commands()->append(&_eventInfoList);
  commands()->append(&_btaTruthList);
  commands()->append(&_applyBoost);
  commands()->append(&_pepCollisionKey);
  commands()->append(&_makeNtuple);
  commands()->append(&_numPrintEvents);
}

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

// The destructor should be limited to undoing the work of the constructor
GqaSemileptonic::~GqaSemileptonic( )
{
}

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

// The beginJob member function is run before any events are
// processed.  In this analysis, it opens the output histogram file
// and then books a number of histograms and a ntuple.

AppResult
GqaSemileptonic::beginJob( AbsEvent* anEvent )
{
  ErrMsg(routine) << name()<< " beginJob "<< endmsg; 
  
  HepTupleManager* manager = gblEnv->getGen()->ntupleManager();
  assert(manager != 0);
  
  _analysisTuple = manager->ntuple("GqaSemileptonic Tuple");

  return AppResult::OK;
}


//--------------------------------------------------------------------
AppResult
GqaSemileptonic::endJob( AbsEvent* anEvent )
{
  ErrMsg(routine) << name( ) << " endJob" << endmsg;
  return AppResult::OK;
}
 


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

AppResult
GqaSemileptonic::event(AbsEvent* anEvent){
  ++_evtNum;   // count events
  
  if ( !_makeNtuple.value() ) return AppResult::OK;
  if ((_evtNum%1000)==0) ErrMsg(routine) << "Processing event:"
					 << _evtNum << endmsg;  

  ErrMsg(trace) << "start GqaSemileptonic " << endmsg;
  //  static const double THMIN_CEN = 0.46;
  //  static const double THMAX_CEN = 2.10;

  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 != 0) nTruth = truthList->length();
  ErrMsg(trace) << "nTruth = " << nTruth << endmsg;


  //HepLorentzVector p4CM = eventInfo->cmFrame();
  //_booster.setLorentzVector(p4CM);
  
  Hep3Vector away = -gblEnv->getPep()->pepBeams()->boostCMtoLab();

  // ###########################################################
  // variables to go into ntuples
  int B1ID(0);
  int B1DecMod(0);
  int L1ID(0); // lepton
  int N1ID(0); // neutrino
  int D1ID(0); // direct D dughter of B, can be D, D* or D**
  int P1ID(0); // non resonant pi
  int D1DecMod(0);
  int finalD1ID(0); // final D daughter of B, not D* or D**
  int finalD1DecMod(0);
  double B1px(0), B1py(0), B1pz(0), B1E(0); //CM frame
  double B1p(0); //CM frame
  double B1pLAB(0); //LABframe
  double L1px(0), L1py(0), L1pz(0), L1E(0); //CM frame
  double L1p(0); //CM frame
  double L1pLAB(0); //LABframe
  double D1px(0), D1py(0), D1pz(0), D1E(0); //CM frame
  double D1p(0); //CM frame
  double D1pLAB(0); //LABframe
  double finalD1px(0), finalD1py(0), finalD1pz(0), finalD1E(0); //CM frame
  double finalD1p(0); //CM frame
  double finalD1pLAB(0); //LABframe
  double mcW1(0);
  double cosL1(-2), cosV1(-2), cosChi1(-2);
  double cosBY1(10);

  int B2ID(0);
  int B2DecMod(0);
  int L2ID(0);
  int N2ID(0);
  int D2ID(0);
  int P2ID(0);
  int D2DecMod(0);
  int finalD2ID(0);
  int finalD2DecMod(0);
  double B2px(0), B2py(0), B2pz(0), B2E(0); //CM frame
  double B2p(0); //CM frame
  double B2pLAB(0); //LABframe
  double L2px(0), L2py(0), L2pz(0), L2E(0); //CM frame
  double L2p(0); //CM frame
  double L2pLAB(0); //LABframe
  double D2px(0), D2py(0), D2pz(0), D2E(0); //CM frame
  double D2p(0); //CM frame
  double D2pLAB(0); //LABframe
  double finalD2px(0), finalD2py(0), finalD2pz(0), finalD2E(0); //CM frame
  double finalD2p(0); //CM frame
  double finalD2pLAB(0); //LABframe
  double mcW2(0);
  double cosL2(-2), cosV2(-2), cosChi2(-2);
  double cosBY2(10);

  // variables for decay modes
  int dig1(0), dig2(0), dig3(0), dig4(0), dig5(0), dig6(0), dig7(0);

  
  // loop over truth list to find the B parents
  HepAListIterator<BtaCandidate> iter(*(truthList));
  BtaCandidate* mcB1(0);
  BtaCandidate* mcB2(0);
  BtaCandidate* cand(0);
  while(cand = iter()){
    if(
       cand->pdtEntry()->lundId()!=PdtLund::B_plus &&
       cand->pdtEntry()->lundId()!=PdtLund::B_minus &&
       cand->pdtEntry()->lundId()!=PdtLund::B0 &&
       cand->pdtEntry()->lundId()!=PdtLund::anti_B0) continue;
    if(mcB1==0) {
      mcB1 = cand;
    } else if(mcB2==0) {
      mcB2 = cand;
    } else break;
  }
  if ( mcB1==0 || mcB2==0 ) return AppResult::OK;
  ErrMsg(trace) << "Found B candidate(s) " << mcB1 << " " << mcB2 << endmsg;

  // loop over the daughter of the mcB1
  dig1=0; dig2=0; dig3=0; dig4=0; dig5=0; dig6=0; dig7=0;
  if (mcB1!=0) {
    B1ID = mcB1->pdtEntry()->lundId();
    ErrMsg(trace) << "B1ID = " << B1ID << endmsg;
    if (B1ID == PdtLund::B_plus) dig1 = 1;
    if (B1ID == PdtLund::B_minus) dig1 = 2;
    if (B1ID == PdtLund::B0) dig1 = 3;
    if (B1ID == PdtLund::anti_B0) dig1 = 4;
    HepAListIterator<BtaCandidate> iterDaut1(mcB1->daughterIterator());
    BtaCandidate* daut1(0);
    BtaCandidate* theL1(0);
    BtaCandidate* theD1(0);
    BtaCandidate* theDdautD1(0);
    BtaCandidate* theFinalD1(0);
    int foundD1(0), foundL1(0), foundN1(0);
    ErrMsg(trace) << "looping over B1 daughters" << endmsg;
    while(daut1 = iterDaut1()) {
      int daughtID = daut1->pdtEntry()->lundId();
      if (daughtID == PdtLund::D_plus) 
	{theD1=daut1; theFinalD1=daut1; D1ID=daughtID; foundD1=1; dig3=1;}
      if (daughtID == PdtLund::D_minus) 
	{theD1=daut1; theFinalD1=daut1; D1ID=daughtID; foundD1=1; dig3=2;}
      if (daughtID == PdtLund::D0) 
	{theD1=daut1; theFinalD1=daut1; D1ID=daughtID; foundD1=1; dig3=3;}
      if (daughtID == PdtLund::anti_D0) 
	{theD1=daut1; theFinalD1=daut1; D1ID=daughtID; foundD1=1; dig3=4;}

      if (daughtID == PdtLund::D_star_plus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=1; dig3=1;}
      if (daughtID == PdtLund::D_star_minus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=1; dig3=2;}
      if (daughtID == PdtLund::D_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=1; dig3=3;}
      if (daughtID == PdtLund::anti_D_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=1; dig3=4;}

      if (daughtID == PdtLund::D_0_star_plus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=2; dig3=1;}
      if (daughtID == PdtLund::D_0_star_minus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=2; dig3=2;}
      if (daughtID == PdtLund::D_0_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=2; dig3=3;}
      if (daughtID == PdtLund::anti_D_0_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=2; dig3=4;}

      if (daughtID == PdtLund::D_1_plus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=3; dig3=1;}
      if (daughtID == PdtLund::D_1_minus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=3; dig3=2;}
      if (daughtID == PdtLund::D_10) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=3; dig3=3;}
      if (daughtID == PdtLund::anti_D_10) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=3; dig3=4;}

      if (daughtID == PdtLund::D_2_star_plus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=4; dig3=1;}
      if (daughtID == PdtLund::D_2_star_minus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=4; dig3=2;}
      if (daughtID == PdtLund::D_2_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=4; dig3=3;}
      if (daughtID == PdtLund::anti_D_2_star0) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=4; dig3=4;}

      if (daughtID == PdtLund::D_prime_1_plus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=5; dig3=1;}
      if (daughtID == PdtLund::D_prime_1_minus) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=5; dig3=2;}
      if (daughtID == PdtLund::D_prime_10) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=5; dig3=3;}
      if (daughtID == PdtLund::anti_D_prime_10) 
	{theD1=daut1; D1ID=daughtID; foundD1=1; dig2=5; dig3=4;}

      if (daughtID == PdtLund::e_minus) 
	{L1ID=daughtID; theL1=daut1; foundL1=1; dig4=1;}
      if (daughtID == PdtLund::e_plus) 
	{L1ID=daughtID; theL1=daut1; foundL1=1; dig4=2;}
      if (daughtID == PdtLund::mu_minus) 
	{L1ID=daughtID; theL1=daut1; foundL1=1; dig4=3;}
      if (daughtID == PdtLund::mu_plus)	
	{L1ID=daughtID; theL1=daut1; foundL1=1; dig4=4;}

      if (daughtID == PdtLund::nu_e) {N1ID=daughtID; foundN1=1; dig5=1;}
      if (daughtID == PdtLund::anti_nu_e) {N1ID=daughtID; foundN1=1; dig5=2;}
      if (daughtID == PdtLund::nu_mu) {N1ID=daughtID; foundN1=1; dig5=3;}
      if (daughtID == PdtLund::anti_nu_mu) {N1ID=daughtID; foundN1=1; dig5=4;}

      if (daughtID == PdtLund::pi0) {P1ID=daughtID; dig6=1;}
      if (daughtID == PdtLund::pi_plus) {P1ID=daughtID; dig6=2;}
      if (daughtID == PdtLund::pi_minus) {P1ID=daughtID; dig6=3;}
    }
    if (foundD1 && foundL1 && foundN1) {
      B1DecMod 
	= dig1*100000 + dig2*10000 + dig3*1000 + dig4*100 + dig5*10 + dig6;
    }
    ErrMsg(trace) << "L1ID = " << L1ID << endmsg;
    ErrMsg(trace) << "N1ID = " << N1ID << endmsg;
    ErrMsg(trace) << "D1ID = " << D1ID << endmsg;
    ErrMsg(trace) << "P1ID = " << P1ID << endmsg;
    ErrMsg(trace) << "B1DecMod = " << B1DecMod << endmsg;
    
    // loop over the daughter of the D1 to get D1 decay mode
    if (dig2>0) {
      dig1=dig2; dig2=dig3; 
      dig3=0; dig4=0; dig5=0; dig6=0; dig7=0;
      HepAListIterator<BtaCandidate> iterD1Daut(theD1->daughterIterator());
      BtaCandidate* D1daut(0);
      BtaCandidate* theD1daut(0);
      ErrMsg(trace) << "looping over D1 daughters" << endmsg;
      while(D1daut = iterD1Daut()) {
	int DdaughtID = D1daut->pdtEntry()->lundId();
	if (DdaughtID == PdtLund::D_plus)  
	  {theDdautD1=D1daut; theFinalD1=D1daut; dig3=1;}
	if (DdaughtID == PdtLund::D_minus) 
	  {theDdautD1=D1daut; theFinalD1=D1daut; dig3=2;}
	if (DdaughtID == PdtLund::D0)      
	  {theDdautD1=D1daut; theFinalD1=D1daut; dig3=3;}
	if (DdaughtID == PdtLund::anti_D0) 
	  {theDdautD1=D1daut; theFinalD1=D1daut; dig3=4;}

	if (DdaughtID == PdtLund::D_star_plus)  
	  {theDdautD1=D1daut; theD1daut=D1daut; dig3=5;}
	if (DdaughtID == PdtLund::D_star_minus) 
	  {theDdautD1=D1daut; theD1daut=D1daut; dig3=6;}
	if (DdaughtID == PdtLund::D_star0)      
	  {theDdautD1=D1daut; theD1daut=D1daut; dig3=7;}
	if (DdaughtID == PdtLund::anti_D_star0) 
	  {theDdautD1=D1daut; theD1daut=D1daut; dig3=8;}
	
	if (DdaughtID == PdtLund::pi0)      {dig4=1;}
	if (DdaughtID == PdtLund::pi_plus)  {dig4=2;}
	if (DdaughtID == PdtLund::pi_minus) {dig4=3;}
	if (DdaughtID == PdtLund::gamma)    {dig4=4;}
      } 
    
      // loop over the daughter of the D1daught
      if (dig3>4) {
	HepAListIterator<BtaCandidate> 
	  iterFD1Daut(theD1daut->daughterIterator());
	BtaCandidate* FD1daut(0);
	while(FD1daut = iterFD1Daut()) {
	  int FDdaughtID = FD1daut->pdtEntry()->lundId();
	  if (FDdaughtID == PdtLund::D_plus)  {theFinalD1=FD1daut; dig5=1;}
	  if (FDdaughtID == PdtLund::D_minus) {theFinalD1=FD1daut; dig5=2;}
	  if (FDdaughtID == PdtLund::D0)      {theFinalD1=FD1daut; dig5=3;}
	  if (FDdaughtID == PdtLund::anti_D0) {theFinalD1=FD1daut; dig5=4;}

	  if (FDdaughtID == PdtLund::pi0)     {dig6=1;}
	  if (FDdaughtID == PdtLund::pi_plus) {dig6=2;}
	  if (FDdaughtID == PdtLund::pi_minus){dig6=3;}
	  if (FDdaughtID == PdtLund::gamma)   {dig6=4;}
	}
      }
      D1DecMod
	= dig1*100000 + dig2*10000 + dig3*1000 + dig4*100 + dig5*10 + dig6;
    }
    ErrMsg(trace) << "D1DecMod = " << D1DecMod << endmsg;

    
    // loop over the daughter of the finalD1 to get finalD decay mode
    dig1=0; dig2=0; dig3=0; dig4=0; dig5=0; dig6=0; dig7=0;

    finalD1ID = theFinalD1->pdtEntry()->lundId();
    if (finalD1ID == PdtLund::D_plus)  {dig1=1;}
    if (finalD1ID == PdtLund::D_minus) {dig1=2;}
    if (finalD1ID == PdtLund::D0)      {dig1=3;}
    if (finalD1ID == PdtLund::anti_D0) {dig1=4;}

    HepAListIterator<BtaCandidate>  
      iterFinalD1(theFinalD1->daughterIterator());
    BtaCandidate* finalD1daut(0);
    while(finalD1daut = iterFinalD1()) {
      int finalD1dautID = finalD1daut->pdtEntry()->lundId();
      if (finalD1dautID == PdtLund::K0)      {dig2=1;}
      if (finalD1dautID == PdtLund::K_S0)    {dig2=1;}
      if (finalD1dautID == PdtLund::K_L0)    {dig2=1;}
      if (finalD1dautID == PdtLund::K_plus)  {dig2=2;}
      if (finalD1dautID == PdtLund::K_minus) {dig2=3;}
      if (finalD1dautID == PdtLund::pi0) {
	if (dig3 == 0 ) {
	  dig3 = 1;
	} else if (dig4 == 0) {
	  dig4 = 1;
	} else {
	  dig5 = 1;
	}
      }
      if (finalD1dautID == PdtLund::pi_plus) {
	if (dig3 == 0 ) {
	  dig3 = 2;
	} else if (dig4 == 0) {
	  dig4 = 2;
	} else {
	  dig5 = 2;
	}
      }
      if (finalD1dautID == PdtLund::pi_minus) {
	if (dig3 == 0 ) {
	  dig3 = 3;
	} else if (dig4 == 0) {
	  dig4 = 3;
	} else {
	  dig5 = 3;
	}
      }
    }
    finalD1DecMod
	= dig1*10000 + dig2*1000 + dig3*100 + dig4*10 + dig5;
    ErrMsg(trace) << "finalD1DecMod = " << finalD1DecMod << endmsg;

    // Get momentum
    HepLorentzVector p4B1 = mcB1->p4(); 
    B1pLAB = p4B1.rho();
    p4B1.boost(away); //CM frame 
    B1px = p4B1.x();
    B1py = p4B1.y();
    B1pz = p4B1.z();
    B1E = p4B1.e();
    B1p = p4B1.rho();
    ErrMsg(trace) << "B1px = " << B1px << endmsg;
    ErrMsg(trace) << "B1py = " << B1py << endmsg;
    ErrMsg(trace) << "B1pz = " << B1pz << endmsg;
    ErrMsg(trace) << "B1E = " << B1E << endmsg;
    ErrMsg(trace) << "B1p = " << B1p << endmsg;
    ErrMsg(trace) << "B1pLAB = " << B1pLAB << endmsg;

    HepLorentzVector p4L1 = theL1->p4(); 
    L1pLAB = p4L1.rho();
    p4L1.boost(away); //CM frame 
    L1px = p4L1.x();
    L1py = p4L1.y();
    L1pz = p4L1.z();
    L1E = p4L1.e();
    L1p = p4L1.rho();

    HepLorentzVector p4D1 = theD1->p4(); 
    D1pLAB = p4D1.rho();
    p4D1.boost(away); //CM frame 
    D1px = p4D1.x();
    D1py = p4D1.y();
    D1pz = p4D1.z();
    D1E = p4D1.e();
    D1p = p4D1.rho();

    HepLorentzVector p4finalD1 = theFinalD1->p4(); 
    finalD1pLAB = p4finalD1.rho();
    p4finalD1.boost(away); //CM frame 
    finalD1px = p4finalD1.x();
    finalD1py = p4finalD1.y();
    finalD1pz = p4finalD1.z();
    finalD1E = p4finalD1.e();
    finalD1p = p4finalD1.rho();


     // Calculate w
    double B1m = p4B1.m();
    double D1m = p4D1.m();
    if (B1m>0 && D1m>0) {
      mcW1 = (B1E*D1E - B1px*D1px - B1py*D1py - B1pz*D1pz) / (B1m*D1m);
    }
    ErrMsg(trace) << "mcW1 = " << mcW1 << endmsg;

    // Calculate angles for D*lnu FF
    HepLorentzVector BP4CMS = p4B1;
    HepLorentzVector leptonP4CMS = p4L1;
    HepLorentzVector DstarP4CMS = p4D1;
    HepLorentzVector WP4CMS;
    WP4CMS.setVect(-Hep3Vector(DstarP4CMS));
    WP4CMS.setE(B1m - DstarP4CMS.e());
    
    HepLorentzVector DstarP4BFrame = DstarP4CMS;
    DstarP4BFrame.boost(-BP4CMS.boostVector());
    HepLorentzVector WP4BFrame = WP4CMS;
    WP4BFrame.boost(-BP4CMS.boostVector());
    HepLorentzVector leptonP4WFrame = leptonP4CMS;
    if (WP4CMS.e() >= WP4CMS.rho() ) {
      leptonP4WFrame.boost(-WP4CMS.boostVector());
      
      double Wp = WP4BFrame.rho();
      double Wpx = WP4BFrame.px();
      double Wpy = WP4BFrame.py();
      double Wpz = WP4BFrame.pz();
      double Dstarp = DstarP4BFrame.rho();
      double Dstarpx = DstarP4BFrame.px();
      double Dstarpy = DstarP4BFrame.py();
      double Dstarpz = DstarP4BFrame.pz();
      double Lp = leptonP4WFrame.rho();
      double Lpx = leptonP4WFrame.px();
      double Lpy = leptonP4WFrame.py();
      double Lpz = leptonP4WFrame.pz();
      
      double cosLNum = Wpx*Lpx + Wpy*Lpy + Wpz*Lpz;
      double cosLDen = Wp*Lp;
      if (cosLDen!=0) {
	cosL1 = cosLNum / cosLDen;
      }

      if (theDdautD1 != 0) {
	HepLorentzVector p4DdautD1 = theDdautD1->p4();
	p4DdautD1.boost(away); //CM frame 
	
	HepLorentzVector DP4CMS = p4DdautD1;
	HepLorentzVector DP4DstarFrame = DP4CMS;
	DP4DstarFrame.boost(-DstarP4CMS.boostVector());
	double DautDp = DP4DstarFrame.rho();
	double DautDpx = DP4DstarFrame.px();
	double DautDpy = DP4DstarFrame.py();
	double DautDpz = DP4DstarFrame.pz();

	double cosVNum 
	  = Dstarpx*DautDpx + Dstarpy*DautDpy + Dstarpz*DautDpz;
	double cosVDen = Dstarp*DautDp;
	if (cosVDen!=0) {
	  cosV1 = cosVNum / cosVDen;
	}
	
	double modLpx = Lpx - Wpx * Lp * cosL1 / Wp;
	double modLpy = Lpy - Wpy * Lp * cosL1 / Wp;
	double modLpz = Lpz - Wpz * Lp * cosL1 / Wp;
	double modDautDpx 
	  = DautDpx - Dstarpx * DautDp * cosV1 / Dstarp;
	double modDautDpy 
	  = DautDpy - Dstarpy * DautDp * cosV1 / Dstarp;
	double modDautDpz 
	  = DautDpz - Dstarpz * DautDp * cosV1 / Dstarp;
	double cosChiNum 
	  = modLpx*modDautDpx + modLpy*modDautDpy
	  + modLpz*modDautDpz;
	double cosChiDen
	  = sqrt(modLpx*modLpx + modLpy*modLpy + modLpz*modLpz)*
	  sqrt(modDautDpx*modDautDpx + modDautDpy*modDautDpy
	       + modDautDpz*modDautDpz);
	if (cosChiDen!=0) {
	  cosChi1 = cosChiNum / cosChiDen;
	}
      }
    }
    ErrMsg(trace) << "cosL1 = " << cosL1 << endmsg;
    ErrMsg(trace) << "cosV1 = " << cosV1 << endmsg;
    ErrMsg(trace) << "cosChi1 = " << cosChi1 << endmsg;
    
    // Calculate cosBY
    HepLorentzVector p4Dl1 = p4finalD1 + p4L1;
    double Dl1E = p4Dl1.e();
    double Dl1p = p4Dl1.rho();
    if (B1p!=0 && Dl1p!=0) {
      cosBY1 = (2*B1E*Dl1E - p4B1.m2() - p4Dl1.m2()) / (2*B1p*Dl1p);
    }
    ErrMsg(trace) << "cosBY1 = " << cosBY1 << endmsg;

  }    


  if (mcB2!=0) {
    B2ID = mcB2->pdtEntry()->lundId();
    ErrMsg(trace) << "B2ID = " << B2ID << endmsg;
    if (B2ID == PdtLund::B_plus) dig1 = 1;
    if (B2ID == PdtLund::B_minus) dig1 = 2;
    if (B2ID == PdtLund::B0) dig1 = 3;
    if (B2ID == PdtLund::anti_B0) dig1 = 4;
    HepAListIterator<BtaCandidate> iterDaut2(mcB2->daughterIterator());
    BtaCandidate* daut2(0);
    BtaCandidate* theL2(0);
    BtaCandidate* theD2(0);
    BtaCandidate* theDdautD2(0);
    BtaCandidate* theFinalD2(0);
    int foundD2(0), foundL2(0), foundN2(0);
    ErrMsg(trace) << "looping over B2 daughters" << endmsg;
    while(daut2 = iterDaut2()) {
      int daughtID = daut2->pdtEntry()->lundId();
      if (daughtID == PdtLund::D_plus) 
	{theD2=daut2; theFinalD2=daut2; D2ID=daughtID; foundD2=1; dig3=1;}
      if (daughtID == PdtLund::D_minus) 
	{theD2=daut2; theFinalD2=daut2; D2ID=daughtID; foundD2=1; dig3=2;}
      if (daughtID == PdtLund::D0) 
	{theD2=daut2; theFinalD2=daut2; D2ID=daughtID; foundD2=1; dig3=3;}
      if (daughtID == PdtLund::anti_D0) 
	{theD2=daut2; theFinalD2=daut2; D2ID=daughtID; foundD2=1; dig3=4;}

      if (daughtID == PdtLund::D_star_plus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=1; dig3=1;}
      if (daughtID == PdtLund::D_star_minus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=1; dig3=2;}
      if (daughtID == PdtLund::D_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=1; dig3=3;}
      if (daughtID == PdtLund::anti_D_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=1; dig3=4;}

      if (daughtID == PdtLund::D_0_star_plus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=2; dig3=1;}
      if (daughtID == PdtLund::D_0_star_minus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=2; dig3=2;}
      if (daughtID == PdtLund::D_0_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=2; dig3=3;}
      if (daughtID == PdtLund::anti_D_0_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=2; dig3=4;}

      if (daughtID == PdtLund::D_1_plus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=3; dig3=1;}
      if (daughtID == PdtLund::D_1_minus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=3; dig3=2;}
      if (daughtID == PdtLund::D_10) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=3; dig3=3;}
      if (daughtID == PdtLund::anti_D_10) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=3; dig3=4;}

      if (daughtID == PdtLund::D_2_star_plus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=4; dig3=1;}
      if (daughtID == PdtLund::D_2_star_minus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=4; dig3=2;}
      if (daughtID == PdtLund::D_2_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=4; dig3=3;}
      if (daughtID == PdtLund::anti_D_2_star0) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=4; dig3=4;}

      if (daughtID == PdtLund::D_prime_1_plus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=5; dig3=1;}
      if (daughtID == PdtLund::D_prime_1_minus) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=5; dig3=2;}
      if (daughtID == PdtLund::D_prime_10) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=5; dig3=3;}
      if (daughtID == PdtLund::anti_D_prime_10) 
	{theD2=daut2; D2ID=daughtID; foundD2=1; dig2=5; dig3=4;}

      if (daughtID == PdtLund::e_minus) 
	{L2ID=daughtID; theL2=daut2; foundL2=1; dig4=1;}
      if (daughtID == PdtLund::e_plus) 
	{L2ID=daughtID; theL2=daut2; foundL2=1; dig4=2;}
      if (daughtID == PdtLund::mu_minus) 
	{L2ID=daughtID; theL2=daut2; foundL2=1; dig4=3;}
      if (daughtID == PdtLund::mu_plus)	
	{L2ID=daughtID; theL2=daut2; foundL2=1; dig4=4;}

      if (daughtID == PdtLund::nu_e) {N2ID=daughtID; foundN2=1; dig5=1;}
      if (daughtID == PdtLund::anti_nu_e) {N2ID=daughtID; foundN2=1; dig5=2;}
      if (daughtID == PdtLund::nu_mu) {N2ID=daughtID; foundN2=1; dig5=3;}
      if (daughtID == PdtLund::anti_nu_mu) {N2ID=daughtID; foundN2=1; dig5=4;}

      if (daughtID == PdtLund::pi0) {P2ID=daughtID; dig6=1;}
      if (daughtID == PdtLund::pi_plus) {P2ID=daughtID; dig6=2;}
      if (daughtID == PdtLund::pi_minus) {P2ID=daughtID; dig6=3;}
    }
    if (foundD2 && foundL2 && foundN2) {
      B2DecMod 
	= dig1*100000 + dig2*10000 + dig3*1000 + dig4*100 + dig5*10 + dig6;
    }
    ErrMsg(trace) << "L2ID = " << L2ID << endmsg;
    ErrMsg(trace) << "N2ID = " << N2ID << endmsg;
    ErrMsg(trace) << "D2ID = " << D2ID << endmsg;
    ErrMsg(trace) << "P2ID = " << P2ID << endmsg;
    ErrMsg(trace) << "B2DecMod = " << B2DecMod << endmsg;
    
    // loop over the daughter of the D2 to get D2 decay mode
    if (dig2>0) {
      dig1=dig2; dig2=dig3; 
      dig3=0; dig4=0; dig5=0; dig6=0; dig7=0;
      HepAListIterator<BtaCandidate> iterD2Daut(theD2->daughterIterator());
      BtaCandidate* D2daut(0);
      BtaCandidate* theD2daut(0);
      ErrMsg(trace) << "looping over D2 daughters" << endmsg;
      while(D2daut = iterD2Daut()) {
	int DdaughtID = D2daut->pdtEntry()->lundId();
	if (DdaughtID == PdtLund::D_plus)  
	  {theDdautD2=D2daut; theFinalD2=D2daut; dig3=1;}
	if (DdaughtID == PdtLund::D_minus) 
	  {theDdautD2=D2daut; theFinalD2=D2daut; dig3=2;}
	if (DdaughtID == PdtLund::D0)      
	  {theDdautD2=D2daut; theFinalD2=D2daut; dig3=3;}
	if (DdaughtID == PdtLund::anti_D0) 
	  {theDdautD2=D2daut; theFinalD2=D2daut; dig3=4;}

	if (DdaughtID == PdtLund::D_star_plus)  
	  {theDdautD2=D2daut; theD2daut=D2daut; dig3=5;}
	if (DdaughtID == PdtLund::D_star_minus) 
	  {theDdautD2=D2daut; theD2daut=D2daut; dig3=6;}
	if (DdaughtID == PdtLund::D_star0)      
	  {theDdautD2=D2daut; theD2daut=D2daut; dig3=7;}
	if (DdaughtID == PdtLund::anti_D_star0) 
	  {theDdautD2=D2daut; theD2daut=D2daut; dig3=8;}
	
	if (DdaughtID == PdtLund::pi0)      {dig4=1;}
	if (DdaughtID == PdtLund::pi_plus)  {dig4=2;}
	if (DdaughtID == PdtLund::pi_minus) {dig4=3;}
	if (DdaughtID == PdtLund::gamma)    {dig4=4;}
      } 
    
      // loop over the daughter of the D2daught
      if (dig3>4) {
	HepAListIterator<BtaCandidate> 
	  iterFD2Daut(theD2daut->daughterIterator());
	BtaCandidate* FD2daut(0);
	while(FD2daut = iterFD2Daut()) {
	  int FDdaughtID = FD2daut->pdtEntry()->lundId();
	  if (FDdaughtID == PdtLund::D_plus)  {theFinalD2=FD2daut; dig5=1;}
	  if (FDdaughtID == PdtLund::D_minus) {theFinalD2=FD2daut; dig5=2;}
	  if (FDdaughtID == PdtLund::D0)      {theFinalD2=FD2daut; dig5=3;}
	  if (FDdaughtID == PdtLund::anti_D0) {theFinalD2=FD2daut; dig5=4;}

	  if (FDdaughtID == PdtLund::pi0)     {dig6=1;}
	  if (FDdaughtID == PdtLund::pi_plus) {dig6=2;}
	  if (FDdaughtID == PdtLund::pi_minus){dig6=3;}
	  if (FDdaughtID == PdtLund::gamma)   {dig6=4;}
	}
      }
      D2DecMod
	= dig1*100000 + dig2*10000 + dig3*1000 + dig4*100 + dig5*10 + dig6;
    }
    ErrMsg(trace) << "D2DecMod = " << D2DecMod << endmsg;

    
    // loop over the daughter of the finalD2 to get finalD decay mode
    dig1=0; dig2=0; dig3=0; dig4=0; dig5=0; dig6=0; dig7=0;

    finalD2ID = theFinalD2->pdtEntry()->lundId();
    if (finalD2ID == PdtLund::D_plus)  {dig1=1;}
    if (finalD2ID == PdtLund::D_minus) {dig1=2;}
    if (finalD2ID == PdtLund::D0)      {dig1=3;}
    if (finalD2ID == PdtLund::anti_D0) {dig1=4;}

    HepAListIterator<BtaCandidate>  
      iterFinalD2(theFinalD2->daughterIterator());
    BtaCandidate* finalD2daut(0);
    while(finalD2daut = iterFinalD2()) {
      int finalD2dautID = finalD2daut->pdtEntry()->lundId();
      if (finalD2dautID == PdtLund::K0)      {dig2=1;}
      if (finalD2dautID == PdtLund::K_S0)    {dig2=1;}
      if (finalD2dautID == PdtLund::K_L0)    {dig2=1;}
      if (finalD2dautID == PdtLund::K_plus)  {dig2=2;}
      if (finalD2dautID == PdtLund::K_minus) {dig2=3;}
      if (finalD2dautID == PdtLund::pi0) {
	if (dig3 == 0 ) {
	  dig3 = 1;
	} else if (dig4 == 0) {
	  dig4 = 1;
	} else {
	  dig5 = 1;
	}
      }
      if (finalD2dautID == PdtLund::pi_plus) {
	if (dig3 == 0 ) {
	  dig3 = 2;
	} else if (dig4 == 0) {
	  dig4 = 2;
	} else {
	  dig5 = 2;
	}
      }
      if (finalD2dautID == PdtLund::pi_minus) {
	if (dig3 == 0 ) {
	  dig3 = 3;
	} else if (dig4 == 0) {
	  dig4 = 3;
	} else {
	  dig5 = 3;
	}
      }
    }
    finalD2DecMod
	= dig1*10000 + dig2*1000 + dig3*100 + dig4*10 + dig5;
    ErrMsg(trace) << "finalD2DecMod = " << finalD2DecMod << endmsg;

    // Get momentum
    HepLorentzVector p4B2 = mcB2->p4(); 
    B2pLAB = p4B2.rho();
    p4B2.boost(away); //CM frame 
    B2px = p4B2.x();
    B2py = p4B2.y();
    B2pz = p4B2.z();
    B2E = p4B2.e();
    B2p = p4B2.rho();
    ErrMsg(trace) << "B2px = " << B2px << endmsg;
    ErrMsg(trace) << "B2py = " << B2py << endmsg;
    ErrMsg(trace) << "B2pz = " << B2pz << endmsg;
    ErrMsg(trace) << "B2E = " << B2E << endmsg;
    ErrMsg(trace) << "B2p = " << B2p << endmsg;
    ErrMsg(trace) << "B2pLAB = " << B2pLAB << endmsg;

    HepLorentzVector p4L2 = theL2->p4(); 
    L2pLAB = p4L2.rho();
    p4L2.boost(away); //CM frame 
    L2px = p4L2.x();
    L2py = p4L2.y();
    L2pz = p4L2.z();
    L2E = p4L2.e();
    L2p = p4L2.rho();

    HepLorentzVector p4D2 = theD2->p4(); 
    D2pLAB = p4D2.rho();
    p4D2.boost(away); //CM frame 
    D2px = p4D2.x();
    D2py = p4D2.y();
    D2pz = p4D2.z();
    D2E = p4D2.e();
    D2p = p4D2.rho();

    HepLorentzVector p4finalD2 = theFinalD2->p4(); 
    finalD2pLAB = p4finalD2.rho();
    p4finalD2.boost(away); //CM frame 
    finalD2px = p4finalD2.x();
    finalD2py = p4finalD2.y();
    finalD2pz = p4finalD2.z();
    finalD2E = p4finalD2.e();
    finalD2p = p4finalD2.rho();


     // Calculate w
    double B2m = p4B2.m();
    double D2m = p4D2.m();
    if (B2m>0 && D2m>0) {
      mcW2 = (B2E*D2E - B2px*D2px - B2py*D2py - B2pz*D2pz) / (B2m*D2m);
    }
    ErrMsg(trace) << "mcW2 = " << mcW2 << endmsg;

    // Calculate angles for D*lnu FF
    HepLorentzVector BP4CMS = p4B2;
    HepLorentzVector leptonP4CMS = p4L2;
    HepLorentzVector DstarP4CMS = p4D2;
    HepLorentzVector WP4CMS;
    WP4CMS.setVect(-Hep3Vector(DstarP4CMS));
    WP4CMS.setE(B2m - DstarP4CMS.e());
    
    HepLorentzVector DstarP4BFrame = DstarP4CMS;
    DstarP4BFrame.boost(-BP4CMS.boostVector());
    HepLorentzVector WP4BFrame = WP4CMS;
    WP4BFrame.boost(-BP4CMS.boostVector());
    HepLorentzVector leptonP4WFrame = leptonP4CMS;
    if (WP4CMS.e() >= WP4CMS.rho() ) {
      leptonP4WFrame.boost(-WP4CMS.boostVector());
      
      double Wp = WP4BFrame.rho();
      double Wpx = WP4BFrame.px();
      double Wpy = WP4BFrame.py();
      double Wpz = WP4BFrame.pz();
      double Dstarp = DstarP4BFrame.rho();
      double Dstarpx = DstarP4BFrame.px();
      double Dstarpy = DstarP4BFrame.py();
      double Dstarpz = DstarP4BFrame.pz();
      double Lp = leptonP4WFrame.rho();
      double Lpx = leptonP4WFrame.px();
      double Lpy = leptonP4WFrame.py();
      double Lpz = leptonP4WFrame.pz();
      
      double cosLNum = Wpx*Lpx + Wpy*Lpy + Wpz*Lpz;
      double cosLDen = Wp*Lp;
      if (cosLDen!=0) {
	cosL2 = cosLNum / cosLDen;
      }

      if (theDdautD2 != 0) {
	HepLorentzVector p4DdautD2 = theDdautD2->p4();
	p4DdautD2.boost(away); //CM frame 
	
	HepLorentzVector DP4CMS = p4DdautD2;
	HepLorentzVector DP4DstarFrame = DP4CMS;
	DP4DstarFrame.boost(-DstarP4CMS.boostVector());
	double DautDp = DP4DstarFrame.rho();
	double DautDpx = DP4DstarFrame.px();
	double DautDpy = DP4DstarFrame.py();
	double DautDpz = DP4DstarFrame.pz();

	double cosVNum 
	  = Dstarpx*DautDpx + Dstarpy*DautDpy + Dstarpz*DautDpz;
	double cosVDen = Dstarp*DautDp;
	if (cosVDen!=0) {
	  cosV2 = cosVNum / cosVDen;
	}
	
	double modLpx = Lpx - Wpx * Lp * cosL2 / Wp;
	double modLpy = Lpy - Wpy * Lp * cosL2 / Wp;
	double modLpz = Lpz - Wpz * Lp * cosL2 / Wp;
	double modDautDpx 
	  = DautDpx - Dstarpx * DautDp * cosV2 / Dstarp;
	double modDautDpy 
	  = DautDpy - Dstarpy * DautDp * cosV2 / Dstarp;
	double modDautDpz 
	  = DautDpz - Dstarpz * DautDp * cosV2 / Dstarp;
	double cosChiNum 
	  = modLpx*modDautDpx + modLpy*modDautDpy
	  + modLpz*modDautDpz;
	double cosChiDen
	  = sqrt(modLpx*modLpx + modLpy*modLpy + modLpz*modLpz)*
	  sqrt(modDautDpx*modDautDpx + modDautDpy*modDautDpy
	       + modDautDpz*modDautDpz);
	if (cosChiDen!=0) {
	  cosChi2 = cosChiNum / cosChiDen;
	}
      }
    }
    ErrMsg(trace) << "cosL2 = " << cosL2 << endmsg;
    ErrMsg(trace) << "cosV2 = " << cosV2 << endmsg;
    ErrMsg(trace) << "cosChi2 = " << cosChi2 << endmsg;
    
    // Calculate cosBY
    HepLorentzVector p4Dl1 = p4finalD2 + p4L2;
    double Dl1E = p4Dl1.e();
    double Dl1p = p4Dl1.rho();
    if (B2p!=0 && Dl1p!=0) {
      cosBY2 = (2*B2E*Dl1E - p4B2.m2() - p4Dl1.m2()) / (2*B2p*Dl1p);
    }
    ErrMsg(trace) << "cosBY2 = " << cosBY2 << endmsg;

  }    





  
  // ############################################################
  // fill ntuple column
  ErrMsg(trace) << "fill ntuple" << endmsg;
  {
    _analysisTuple->column("nMC", nTruth);

    _analysisTuple->column("B1ID", B1ID);
    _analysisTuple->column("B1DecMod", B1DecMod);
    _analysisTuple->column("Lep1ID", L1ID);
    _analysisTuple->column("Neu1ID", N1ID);
    _analysisTuple->column("D1ID", D1ID);
    _analysisTuple->column("Pi1ID", P1ID);
    _analysisTuple->column("D1DecMod", D1DecMod);
    _analysisTuple->column("finalD1ID", finalD1ID);
    _analysisTuple->column("finalD1DecMod", finalD1DecMod);
    _analysisTuple->column("B1px", B1px);
    _analysisTuple->column("B1py", B1py);
    _analysisTuple->column("B1pz", B1pz);
    _analysisTuple->column("B1E", B1E);
    _analysisTuple->column("B1p", B1p);
    _analysisTuple->column("B1pLAB", B1pLAB);
    _analysisTuple->column("L1px", L1px);
    _analysisTuple->column("L1py", L1py);
    _analysisTuple->column("L1pz", L1pz);
    _analysisTuple->column("L1E", L1E);
    _analysisTuple->column("L1p", L1p);
    _analysisTuple->column("L1pLAB", L1pLAB);
    _analysisTuple->column("D1px", D1px);
    _analysisTuple->column("D1py", D1py);
    _analysisTuple->column("D1pz", D1pz);
    _analysisTuple->column("D1E", D1E);
    _analysisTuple->column("D1p", D1p);
    _analysisTuple->column("D1pLAB", D1pLAB);
    _analysisTuple->column("finalD1px", finalD1px);
    _analysisTuple->column("finalD1py", finalD1py);
    _analysisTuple->column("finalD1pz", finalD1pz);
    _analysisTuple->column("finalD1E", finalD1E);
    _analysisTuple->column("finalD1p", finalD1p);
    _analysisTuple->column("finalD1pLAB", finalD1pLAB);
    _analysisTuple->column("mcW1", mcW1);
    _analysisTuple->column("cosL1", cosL1);
    _analysisTuple->column("cosV1", cosV1);
    _analysisTuple->column("cosChi1", cosChi1);
    _analysisTuple->column("cosBY1", cosBY1);

    _analysisTuple->column("B2ID", B2ID);
    _analysisTuple->column("B2DecMod", B2DecMod);
    _analysisTuple->column("Lep2ID", L2ID);
    _analysisTuple->column("Neu2ID", N2ID);
    _analysisTuple->column("D2ID", D2ID);
    _analysisTuple->column("Pi2ID", P2ID);
    _analysisTuple->column("D2DecMod", D2DecMod);
    _analysisTuple->column("finalD2ID", finalD2ID);
    _analysisTuple->column("finalD2DecMod", finalD2DecMod);
    _analysisTuple->column("B2px", B2px);
    _analysisTuple->column("B2py", B2py);
    _analysisTuple->column("B2pz", B2pz);
    _analysisTuple->column("B2E", B2E);
    _analysisTuple->column("B2p", B2p);
    _analysisTuple->column("B2pLAB", B2pLAB);
    _analysisTuple->column("L2px", L2px);
    _analysisTuple->column("L2py", L2py);
    _analysisTuple->column("L2pz", L2pz);
    _analysisTuple->column("L2E", L2E);
    _analysisTuple->column("L2p", L2p);
    _analysisTuple->column("L2pLAB", L2pLAB);
    _analysisTuple->column("D2px", D2px);
    _analysisTuple->column("D2py", D2py);
    _analysisTuple->column("D2pz", D2pz);
    _analysisTuple->column("D2E", D2E);
    _analysisTuple->column("D2p", D2p);
    _analysisTuple->column("D2pLAB", D2pLAB);
    _analysisTuple->column("finalD2px", finalD2px);
    _analysisTuple->column("finalD2py", finalD2py);
    _analysisTuple->column("finalD2pz", finalD2pz);
    _analysisTuple->column("finalD2E", finalD2E);
    _analysisTuple->column("finalD2p", finalD2p);
    _analysisTuple->column("finalD2pLAB", finalD2pLAB);
    _analysisTuple->column("mcW2", mcW2);
    _analysisTuple->column("cosL2", cosL2);
    _analysisTuple->column("cosV2", cosV2);
    _analysisTuple->column("cosChi2", cosChi2);
    _analysisTuple->column("cosBY2", cosBY2);
  }

  // dump to ntuple
  _analysisTuple->dumpData();
  

  return AppResult::OK;
}
 
//--------------------------------------------------------------------