////////////////////////////////////////////////////////////////////////
//                                                                    // 
//Elisabetta 15/08/2013                                               //
// starting macro                                                     //
//Analysis: p pbar ---> Ds1+ Ds-                                      //
//Ds1 = Ds(2460)+                                                     //
//where: Ds1(2460) ---> D*s+ pi0, D*s+ to Ds+(--->K+ K- pi+) gamma    //
//                                                                    //
//                                                                    //
//                                                                    //
////////////////////////////////////////////////////////////////////////
class RhoCandList;
class RhoCandidate;
class PndAnaPidSelector;
class PndAnaPidCombiner;
class PndAnalysis;
#include <iostream>
using std::cout;
using std::endl;


// *** routine to only keep PID matched candidates in list
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
{
	int removed = 0;
	
	for (int ii=l.GetLength()-1;ii>=0;--ii)
	{
		if ( !(ana->McTruthMatch(l[ii])) )
		{
			l.Remove(l[ii]);
			removed++;
		}
	}
	
	return removed;
}


void tut_anaDs2460Dstarpi0_simple(int nevts=5000)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	
	TString OutFile="output.root";  
					
	// *** the files coming from the simulation
	TString inPidFile = "pid_complete.root";    // this file contains the PndPidCandidates and McTruth
	TString inParFile = "simparams.root";
	
	gStyle->SetOptFit(1011);
	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
	
	FairLogger::GetLogger()->SetLogToFile(kFALSE);
	
	// *** initialization
	FairRunAna* fRun = new FairRunAna();
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
	fRun->SetInputFile(inPidFile);
	
	TString anaParFile = gSystem->Getenv("VMCWORKDIR");
	anaParFile +=  "/macro/params/all.par";
	FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
	parIo1->open(anaParFile.Data(),"in");
        

	FairParRootFileIo* parIO = new FairParRootFileIo();
	parIO->open(inParFile);
	rtdb->setFirstInput(parIO);
	rtdb->setSecondInput(parIo1);
	rtdb->setOutput(parIO);  
	

	fRun->SetOutputFile(OutFile);
	fRun->Init(); 
	
	// *** create an output file for all histograms
	TFile *out = TFile::Open("output_anaDs2460Dstarpi0.root","RECREATE");
		
	// *** create some histograms
	TH1F *gam_en = new TH1F("gam_en","#gamma energy",200,0.,5);
	TH1F *pi0_mom = new TH1F("pi0_mom","#pi0 momentum",200,0.,5);

	/*
	TH1F *Dsp_lpid = new TH1F("Dsp_lpid","Ds+ mass (loose pid)",100,0.,5.);
	TH1F *Dsm_lpid = new TH1F("Dsm_lpid","Ds- mass (loose pid)",100,0.,5.);
	TH1F *D2460_lpid = new TH1F("D2460_lpid","Ds1^{+} mass (loose pid)",100,0.,5.);
	TH1F *pi0_lpid  = new TH1F("pi0_lpid","#pi0 mass (loose pid)",100,0.,0.6);
	TH1F *Dstar_lpid  = new TH1F("Dstar_lpid","Ds*+ mass (loose pid)",100,0.,5.);
	*/

	TH1F *Dsp_best = new TH1F("Dsp_best","Ds+ mass (best pid)",100,0.,5.);
	TH1F *Dsm_best = new TH1F("Dsm_best","Ds- mass (best pid)",100,0.,5.);
	TH1F *D2460_best = new TH1F("D2460_best","Ds1^{+} mass (best pid)",100,0.,5.);
	TH1F *pi0_best  = new TH1F("pi0_best","#pi0 mass (best pid)",100,0.,0.6);
	TH1F *Dstar_best  = new TH1F("Dstar_best","Ds*+ mass (best pid)",100,0.,5.);
	TH1F *pi0_ftm = new TH1F("pi0_ftm","#pi0 mass (full truth match)",100,0.,0.6);
	TH1F *pi0_Mom_ftm = new TH1F("pi0_Mom_ftm","#pi0 mass (full truth match)",100,0.,3.);
	TH1F *Dsp_Mom_ftm = new TH1F("Dsp_Mom_ftm","#pi0 mass (full truth match)",100,0.,3.);
	TH1F *Dsm_Mom_ftm = new TH1F("Dsm_Mom_ftm","#pi0 mass (full truth match)",100,0.,3.);
	TH1F *Dstar_Mom_ftm = new TH1F("Dsp_Mom_ftm","#pi0 mass (full truth match)",100,0.,3.);
	TH1F *D2460_Mom_ftm = new TH1F("Dsp_Mom_ftm","#pi0 mass (full truth match)",100,0.,3.);

	//pi0 plots: all and ftm
	TH1F *pi0_all = new TH1F("pi0_all","#pi0 mass (all)",100,0.,0.6);
	TH1F *pi0_px_diff_ftm  = new TH1F("pi0_px_diff_ftm","#pi0 px ",100,-0.5,0.5);
	TH1F *pi0_py_diff_ftm  = new TH1F("pi0_py_diff_ftm","#pi0 py ",100,-0.5,0.5);
	TH1F *pi0_pz_diff_ftm  = new TH1F("pi0_pz_diff_ftm","#pi0 pz ",100,-0.5,0.5);
	TH1F *pi0_E_diff_ftm  = new TH1F("pi0_E_diff_ftm","#pi0 pz ",100,-0.5,0.5);
	TH1F *pi0_mass_diff_ftm  = new TH1F("pi0_mass_diff_ftm","#pi0 mass resolution ",100,-0.25,0.25);
	TH1F *pi0_px_ftm  = new TH1F("pi0_px_ftm","#pi0 px ",100,-1.5,1.5);
	TH1F *pi0_py_ftm  = new TH1F("pi0_py_ftm","#pi0 py ",100,-1.5,1.5);
	TH1F *pi0_pz_ftm  = new TH1F("pi0_pz_ftm","#pi0 pz ",100,-0.5,3.);
	TH1F *pi0_E_ftm  = new TH1F("pi0_E_ftm","#pi0 E ",100,0,10);
	TH1F *pi0_ftm = new TH1F("pi0_ftm","#pi0 mass (full truth match)",100,0.,0.6);

	//Ds*+ plots: all and ftm
	TH1F *Dstar_all = new TH1F("Dstar_all","Ds*+ mass (all)",100,0.,5.);
	TH1F *Dstar_px_diff_ftm  = new TH1F("Dstar_px_diff_ftm","Ds*+ px ",100,-0.5,0.5);
	TH1F *Dstar_py_diff_ftm  = new TH1F("Dstar_py_diff_ftm","Ds*+ py ",100,-0.5,0.5);
	TH1F *Dstar_pz_diff_ftm  = new TH1F("Dstar_pz_diff_ftm","Ds*+ pz ",100,-0.5,0.5);
	TH1F *Dstar_E_diff_ftm  = new TH1F("Dstar_E_diff_ftm","Ds*+ pz ",100,-0.5,0.5);
	TH1F *Dstar_mass_diff_ftm  = new TH1F("Dstar_mass_diff_ftm","Ds*+ mass resolution ",100,-0.25,0.25);
	TH1F *Dstar_px_ftm  = new TH1F("Dstar_px_ftm","Ds*+ px ",100,-1.5,1.5);
	TH1F *Dstar_py_ftm  = new TH1F("Dstar_py_ftm","Ds*+ py ",100,-1.5,1.5);
	TH1F *Dstar_pz_ftm  = new TH1F("Dstar_pz_ftm","Ds*+ pz ",100,-1,9);
	TH1F *Dstar_E_ftm  = new TH1F("Dstar_E_ftm","Ds*+ E ",100,0,10);
	TH1F *Dstar_ftm = new TH1F("Dstar_ftm","Ds*+ mass (full truth match)",100,0.,5.);


	//Ds plots: all and ftm
	TH1F *Dsp_all  = new TH1F("Dsp_all","Ds+ mass (all)",100,0.,5.);
	TH1F *Dsp_px_ftm  = new TH1F("Dsp_px_ftm","Ds+ px ",100,-1.5,1.5);
	TH1F *Dsp_py_ftm  = new TH1F("Dsp_py_ftm","Ds+ py ",100,-1.5,1.5);
	TH1F *Dsp_pz_ftm  = new TH1F("Dsp_pz_ftm","Ds+ pz ",100,-1,9);
	TH1F *Dsp_E_ftm  = new TH1F("Dsp_E_ftm","Ds+ E ",100,0,10);
	TH1F *Dsp_ftm  = new TH1F("Dsp_ftm","Ds+ mass (full truth match)",100,0.,5.);
	TH1F *Dsp_mass_diff_ftm  = new TH1F("Dsp_mass_diff_ftm","Ds+ E ",100,-0.2,0.2);
	TH1F *Dsp_px_diff_ftm  = new TH1F("Dsp_px_diff_ftm","Ds+ px ",100,-0.5,0.5);
	TH1F *Dsp_py_diff_ftm  = new TH1F("Dsp_py_diff_ftm","Ds+ py ",100,-0.5,0.5);
	TH1F *Dsp_pz_diff_ftm  = new TH1F("Dsp_pz_diff_ftm","Ds+ pz ",100,-0.5,0.5);
	TH1F *Dsp_E_diff_ftm  = new TH1F("Dsp_E_ftm","Ds+ E ",100,-0.5,0.5);

	TH1F *Dsm_all  = new TH1F("Dsm_all","Ds- mass (all)",100,0.,5.);
	TH1F *Dsm_px_ftm  = new TH1F("Dsm_px_ftm","Ds- px ",100,-1.5,1.5);
	TH1F *Dsm_py_ftm  = new TH1F("Dsm_py_ftm","Ds- py ",100,-1.5,1.5);
	TH1F *Dsm_pz_ftm  = new TH1F("Dsm_pz_ftm","Ds- pz ",100,-1,9);
	TH1F *Dsm_E_ftm  = new TH1F("Dsm_E_ftm","Ds- E ",100,0,10);
	TH1F *Dsm_ftm  = new TH1F("Dsm_ftm","Ds- mass (full truth match)",100,0,5.);
	TH1F *Dsm_mass_diff_ftm  = new TH1F("Dsm_mass_diff_ftm","Ds- E ",100,-0.2,0.2);
	TH1F *Dsm_px_diff_ftm  = new TH1F("Dsm_px_diff_ftm","Ds- px ",100,-0.5,0.5);
	TH1F *Dsm_py_diff_ftm  = new TH1F("Dsm_py_diff_ftm","Ds- py ",100,-0.5,0.5);
	TH1F *Dsm_pz_diff_ftm  = new TH1F("Dsm_pz_diff_ftm","Ds- pz ",100,-0.5,0.5);
	TH1F *Dsm_E_diff_ftm  = new TH1F("Dsm_E_ftm","Ds- E ",100,-0.5,0.5);


	
	//Ds2460 plots: all and ftm
	TH1F *D2460_all  = new TH1F("D2460_all","Ds'1^{+} mass (all)",100,0.,5.);
	TH1F *D2460_px_ftm  = new TH1F("D2460_px_ftm","Ds'1+ px ",100,-1.5,1.5);
	TH1F *D2460_py_ftm  = new TH1F("D2460_py_ftm","Ds'1+ py ",100,-1.5,1.5);
	TH1F *D2460_pz_ftm  = new TH1F("D2460_pz_ftm","Ds'1+ pz ",100,-1,9);
	TH1F *D2460_E_ftm  = new TH1F("D2460_E_ftm","Ds'1+ E ",100,0,10);
	TH1F *D2460_mass_diff_ftm  = new TH1F("D2460_mass_diff_ftm","Ds'1+ mass resolution (full truth match)",100,-0.2,0.2);
	TH1F *D2460_ftm  = new TH1F("D2460_ftm","Ds'1+ mass (full truth match)",100,0.,4.);
	TH1F *D2460_px_diff_ftm  = new TH1F("D2460_px_diff_ftm","Ds'1+ px resolution",100,-0.5,0.5);
	TH1F *D2460_py_diff_ftm  = new TH1F("D2460_py_diff_ftm","Ds'1+ py resolution",100,-0.5,0.5);
	TH1F *D2460_pz_diff_ftm  = new TH1F("D2460_pz_diff_ftm","Ds'1+ pz resolution",100,-0.5,0.5);
	TH1F *D2460_E_diff_ftm  = new TH1F("D2460_E_diff_ftm","Ds'1+ E resolution",100,-0.5,0.5);
	TH1F *D2460_diff_mcf  = new TH1F("D2460_diff_mcf","Ds'1+ mass resolution",100,-0.5,0.5);

	//best 

	TH1F *hpxd0_best  = new TH1F("hpxd0_best","D0 px ",100,-1.5,1.5);
	TH1F *hpxd0star_best  = new TH1F("hpxd0star_best","D0* px ",100,-1.5,1.5);
	TH1F *hpxds_best  = new TH1F("hpxds_best","Ds px ",100,-1.5,1.5);
	TH1F *hpyd0_best  = new TH1F("hpyd0_best","D0 py ",100,-1.5,1.5);
	TH1F *hpyd0star_best  = new TH1F("hpyd0star_best","D0* py ",100,-1.5,1.5);
	TH1F *hpyds_best  = new TH1F("hpyds_best","Ds py ",100,-1.5,1.5);
	TH1F *hpzd0_best  = new TH1F("hpzd0_best","D0 pz ",100,-1,9);
	TH1F *hpzd0star_best  = new TH1F("hpzd0star_best","D0* pz ",100,-1,9);
	TH1F *hpzds_best  = new TH1F("hpzds_best","Ds pz ",100,-1,9);
	TH1F *hEds_best  = new TH1F("hEds_best","Ds E ",100,0,10);
	TH1F *hEd0_best  = new TH1F("hEd0_best","D0 E ",100,0,10);
	TH1F *hEd0star_best  = new TH1F("hEd0star_best","D0* E ",100,0,10);
	TH1F *hpxd2460_best  = new TH1F("hpxd2460_best","Ds'1+ px ",100,-1.5,1.5);
	TH1F *hpyd2460_best  = new TH1F("hpyd2460_best","Ds'1+ py ",100,-1.5,1.5);
	TH1F *hpzd2460_best  = new TH1F("hpzd2460_best","Ds'1+ pz ",100,-1,9);
	TH1F *hEd2460_best  = new TH1F("hEd2460_best","Ds'1+ E ",100,0,10);
	
	TH1F *hpmpx_m  = new TH1F("hpmpx_m","#pi- px, daughter of Ds- ",100,-3,3);
	TH1F *hpmpy_m  = new TH1F("hpmpy_m","#pi- py, daughter of Ds- ",100,-3,3);
	TH1F *hpmpz_m = new TH1F("hpmpz_m","#pi- pz, daughter of Ds- ",100,-1,9);
	TH1F *hpmE_m = new TH1F("hpmE_m","#pi- E, daughter of Ds- ",100,0,10);

	TH1F *hkmpx_m  = new TH1F("hkmpx_m","K- px, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkmpy_m  = new TH1F("hkmpy_m","K- py, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkmpz_m = new TH1F("hkmpz_m","K- pz, daughter of Ds- ",100,-1,9);
	TH1F *hkmE_m = new TH1F("hkmE_m","K- E, daughter of Ds- ",100,0,10);
	
	TH1F *hkppx_m  = new TH1F("hkppx_m","K+ px, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkppy_m  = new TH1F("hkppy_m","K+ py, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkppz_m = new TH1F("hkppz_m","K+ pz, daughter of Ds- ",100,-1,9);
	TH1F *hkpE_m = new TH1F("hkpE_m","K+ E, daughter of Ds- ",100,0,10);

	TH1F *h12_m  = new TH1F("h12_m","Ds-: K+ K- invariant mass",100,0.5,3);
	TH1F *h13_m  = new TH1F("h13_m","Ds-: K+ pi+ invariant mass",100,0.5,3);
	TH1F *h23_m  = new TH1F("h23_m","Ds-: K- pi+ invariant mass",100,0.5,3);
	

	//

		
	TH1F *hpmpx_p  = new TH1F("hpmpx_p","#pi- px, daughter of Ds- ",100,-3,3);
	TH1F *hpmpy_p  = new TH1F("hpmpy_p","#pi- py, daughter of Ds- ",100,-3,3);
	TH1F *hpmpz_p = new TH1F("hpmpz_p","#pi- pz, daughter of Ds- ",100,-1,9);
	TH1F *hpmE_p = new TH1F("hpmE_p","#pi- E, daughter of Ds- ",100,0,10);

	TH1F *hkmpx_p  = new TH1F("hkmpx_p","K- px, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkmpy_p  = new TH1F("hkmpy_p","K- py, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkmpz_p = new TH1F("hkmpz_p","K- pz, daughter of Ds- ",100,-1,9);
	TH1F *hkmE_p = new TH1F("hkmE_p","K- E, daughter of Ds- ",100,0,10);
	
	TH1F *hkppx_p  = new TH1F("hkppx_p","K+ px, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkppy_p  = new TH1F("hkppy_p","K+ py, daughter of Ds- ",100,-1.5,1.5);
	TH1F *hkppz_p = new TH1F("hkppz_p","K+ pz, daughter of Ds- ",100,-1,9);
	TH1F *hkpE_p = new TH1F("hkpE_p","K+ E, daughter of Ds- ",100,0,10);

	TH1F *h12_p  = new TH1F("h12_p","Ds+: K+ K- invariant mass",100,0.5,3);
	TH1F *h13_p  = new TH1F("h13_p","Ds+: K+ pi+ invariant mass",100,0.5,3);
	TH1F *h23_p  = new TH1F("h23_p","Ds+: K- pi+ invariant mass",100,0.5,3);
	
	TH1F *h12_m  = new TH1F("h12_m","Ds-: K+ K- invariant mass",100,0.5,3);
	TH1F *h13_m  = new TH1F("h13_m","Ds-: K+ pi- invariant mass",100,0.5,3);
	TH1F *h23_m  = new TH1F("h23_m","Ds-: K- pi- invariant mass",100,0.5,3);

	TH1F *pi0_tpid = new TH1F("pi0_tpid","#pi0 mass (tight pid)",100,0,0.6);
	TH1F *Dstar_tpid = new TH1F("Dstar_tpid","Ds*+ mass (tight pid)",100,0,5.);
	TH1F *Dsm_tpid = new TH1F("Dsm_tpid","Ds- mass (tight pid)",100,0,5.);
	TH1F *Dsp_tpid = new TH1F("Dsp_tpid","Ds+ mass (tight pid)",100,0,5.);
	TH1F *D2460_tpid  = new TH1F("D2460_tpid"," Ds1^{+} mass (tight pid)",100,0,5);
	
	TH1F *D2460_trpid = new TH1F("D2460_trpid","Ds'1^{+} mass (true pid)",100,0,5.);
	TH1F *pi0_trpid = new TH1F("pi0_trpid","#pi0 mass (true pid)",100,0,0.6);
	TH1F *Dstar_trpid = new TH1F("Dstar_trpid","Ds*+ mass (true pid)",100,0,5.);
	TH1F *Dsm_trpid  = new TH1F("Dsm_trpid","Ds- mass (true pid)",100,0,5);
	TH1F *Dsp_trpid  = new TH1F("Dsp_trpid","Ds+ mass (true pid)",100,0,5);
	
	
	TH1F *Dsp_nm = new TH1F("Dsp_nm","Ds+ mass (no truth match)",100,0,5.);
	TH1F *Dsm_nm = new TH1F("Dsm_nm","Ds- mass (no truth match)",100,0,5.);
	TH1F *pi0_nm = new TH1F("pi0_nm","#pi0 mass (no truth match)",100,0,0.6);
	TH1F *Dstar_nm = new TH1F("Dstarstar_nm","Ds*+ mass (no truth match)",100,0,5.);
	TH1F *D2460_nm  = new TH1F("D2460_nm","Ds'1+ mass (no truth match)",100,0,5.);
		

	TH1F *hm_diff  = new TH1F("hm_diff","Ds1+ mass diff to truth (4C)",100,0,10);
	TH1F *hm_diff4  = new TH1F("hm_diff4","Ds*+ mass diff to truth after fit",100,-0.5,0.5);
	TH1F *hm_diff2m  = new TH1F("hm_diff2m","Ds- mass diff to truth after vertex fit",100,-0.5,0.5);
	TH1F *hm_diff2p  = new TH1F("hm_diff2p","Ds+ mass diff to truth after vertex fit",100,-0.5,0.5);
	TH1F *hm_diff3  = new TH1F("hm_diff3","#pi0 mass diff to truth after vertex fit",100,-0.5,0.5);
	TH1F *hm_diff2m_mcf  = new TH1F("hm_diff2m_mcf","Ds- mass diff to truth after mcf",100,-0.5,0.5);
	TH1F *hm_diff2p_mcf  = new TH1F("hm_diff2p_mcf","Ds+ mass diff to truth after mcf",100,-0.5,0.5);
	TH1F *hm_diffstar_mcf  = new TH1F("hm_diffstar_mcf","Ds*+ mass diff to truth after mcf",100,-0.5,0.5);
	TH1F *hm_diffpi0_mcf  = new TH1F("hm_diffpi0_mcf","#pi0 mass diff to truth after mcf",100,-0.5,0.5);
	TH1F *hm_diffds1_mcf  = new TH1F("hm_diffds1_mcf","Ds1+ mass diff to truth after mcf",100,-0.5,0.5);
 
	//fitters

	TH1F *Dsp_vf   = new TH1F("Dsp_vf","Ds+ mass (vertex fit)",100,0,6);
	TH1F *Dsm_vf   = new TH1F("Dsm_vf","Ds- mass (vertex fit)",100,0,6);
	TH1F *D2460_mfc   = new TH1F("D2460_mfc","Ds1+ mass (vertex fit)",100,0,6);
	TH1F *pi0_vf   = new TH1F("pi0_vf","pi0 mass (vertex fit)",100,0,6);
	TH1F *Dstar_vf   = new TH1F("Dstar_vf","Ds*+ mass (vertex fit)",100,0,10);

	TH1F *D2460_4cf  = new TH1F("D2460_4cf","Ds1+ mass (4C fit)",100,0,6);
	TH1F *Dsm_4cf  = new TH1F("Dsm_4cf","Ds- mass (4C fit)",100,0,3);
	TH1F *Dsp_4cf  = new TH1F("Dsp_4cf","Ds- mass (4C fit)",100,0,3);
	TH1F *pi0_4cf  = new TH1F("pi0_4cf","#pi0 mass (4C fit)",100,0,0.6);
	TH1F *pi0_mcf  = new TH1F("pi0_mcf","#pi0 mass (mass constraint fit)",100,0,0.6);
	TH1F *Dstar_mcf  = new TH1F("Dstar_mcf","Ds* mass (mass constraint fit)",100,0,5);
	TH1F *Dsp_mcf  = new TH1F("Dsp_mcf","Ds+ mass (mass constraint fit)",100,0,3);
	TH1F *Dsm_mcf  = new TH1F("Dsm_mcf","Ds- mass (mass constraint fit)",100,0,3);
	TH1F *D2460_mcf  = new TH1F("D2460_mcf","D's1+ mass (mass constraint fit)",100,0,6);
	TH1F *Dsp_chi2_vf  = new TH1F("Dsp_chi2_vf", "Ds+: #chi^{2} vertex fit",100,-1,10);
	TH1F *Dsm_chi2_vf  = new TH1F("DsM_chi2_vf", "Ds-: #chi^{2} vertex fit",100,-1,10);
	TH1F *D2460_chi2_mf  = new TH1F("D2460_chi2_mf", "Ds'1: #chi^{2} mass fit",100,-1,10);
	TH1F *pi0_chi2_vf  = new TH1F("pi0_chi2_vf", "#pi0: #chi^{2} vertex fit",100,-1,10);
	TH1F *Dstar_chi2_vf  = new TH1F("Dstar_chi2_vf", "Ds*+: #chi^{2} vertex fit",100,-1,10);
	TH1F *D2460_chi2_vf  = new TH1F("D2460_chi2_vf", "D2460: #chi^{2} vertex fit",100,-1,10);


	TH1F *Dsm_vtx_z   = new TH1F("Dsm_vtx_z",  "Ds- Z",100,-0.5,0.5);
	TH1F *Dsp_vtx_z   = new TH1F("Dsp_vtx_z",  "Ds+ Z",100,-0.5,0.5);
	TH1F *pi0_vtx_z   = new TH1F("pi0_vtx_z",  "#pi0 Z",100,-0.5,0.5);
	TH1F *Dstar_vtx_z   = new TH1F("Dstar_vtx_z",  "Ds*+ Z",100,-0.5,0.5);
	TH1F *D2460_vtx_z   = new TH1F("D2460_vtx_z",  "D'1s+: Z",100,-0.5,0.5);
	TH1F *Dsp_chi2_prob_vf   = new TH1F("Dsp_chi2_prob_vf",  "Ds+ prob vertex fit",100,-0.01,1);
	TH1F *Dsm_chi2_prob_vf   = new TH1F("Dsm_chi2_prob_vf",  "Ds- prob vertex fit",100,-0.01,1);
	TH1F *pi0_chi2_prob_vf   = new TH1F("pi0_chi2_prob_vf",  "#pi0: prob vertex fit",100,-0.01,1);
	TH1F *Dstar_chi2_prob_vf   = new TH1F("Dstar_chi2_prob_vf",  "Ds*: prob vertex fit",100,-0.01,1);
	TH1F *D2460_chi2_prob_vf   = new TH1F("D2460_chi2_prob_vf",  "D's1+: prob vertex fit",100,-0.01,1);
	TH1F *D2460_chi2_prob_mf   = new TH1F("D2460_chi2_prob_mf",  "Ds'1+: prob mass fit",100,-0.01,1);
	TH1F *D2460_chi2_4c   = new TH1F("D2460_chi2_4c",  "Ds'1+: #chi^{2} 4C fit",100,-1,20);
	TH1F *D2460_chi2_prob_4c   = new TH1F("D2460_chi2_prob_4c",  "Ds'1+: prob 4C fit",100,-0.01,1);
	TH1F *pi0_chi2_mf  = new TH1F("pi0_chi2_mf", "#pi0: #chi^{2} mass fit",100,-1,10);
	TH1F *pi0_chi2_prob_mf  = new TH1F("pi0_chi2_prob__mf", "#pi0: prob mass constaint fit",100,0,1);
	TH1F *Dstar_chi2_mf  = new TH1F("Dstar_chi2_mf", "Ds*+: #chi^{2} mass fit",100,-1,10);
	TH1F *Dstar_chi2_vf = new TH1F("Dstar_chi2_vf", "Ds*+: #chi^{2} mass fit",100,-1,10);
	TH1F *Dstar_chi2_prob_mf  = new TH1F("Dstar_chi2_prob_mf", "Ds*+: prob mass constraint fit",100,0,1);

	TH1F *Dsp_chi2_mf  = new TH1F("Dsp_chi2_mf", "Ds-: #chi^{2} mass fit",100,-1,10);
	TH1F *Dsm_chi2_mf  = new TH1F("Dsm_chi2_mf", "Ds+: #chi^{2} mass fit",100,-1,10);
	TH1F *Dsp_chi2_prob_mf  = new TH1F("Dsp_chi2_prob_mf", "Ds+: #chi^{2} mass constrained fit",100,0,1);
	TH1F *Dsm_chi2_prob_mf  = new TH1F("Dsm_chi2_prob_mf", "Ds-: #chi^{2} mass constrained fit",100,0,1);
	
	TH2F *hDalitzm12m13_m = new TH2F("hDalitzm12m13_m","m12 vs m13",30,0,8.5,30,0,8.5);
	TH2F *hDalitzm12m23_m = new TH2F("hDalitzm12m23_m","m12 vs m23",30,0,8.5,30,0,8.5);
	TH2F *hDalitzm13m23_m = new TH2F("hDalitzm13m23_m","m13 vs m23",30,0,8.5,30,0,8.5);
	TH2F *hDalitzm12m13_p = new TH2F("hDalitzm12m13_p","m12 vs m13",30,0,8.5,30,0,8.5);
	TH2F *hDalitzm12m23_p = new TH2F("hDalitzm12m23_p","m12 vs m23",30,0,8.5,30,0,8.5);
	TH2F *hDalitzm13m23_p = new TH2F("hDalitzm13m23_p","m13 vs m23",30,0,8.5,30,0,8.5);
	

	TH2F *hvpos = new TH2F("hvpos","(x,y) projection of ftm Ds*1+",100,-2,2,100,-2,2);
	TH2F *hvpos2p = new TH2F("hvpos2p","(x,y) projection of fitted decay vertex Ds-",100,-2,2,100,-2,2);
	TH2F *hvpos2m = new TH2F("hvpos2m","(x,y) projection of fitted decay vertex Ds-",100,-2,2,100,-2,2);
	TH2F *hvpos3 = new TH2F("hvpos3","(x,y) projection of fitted decay vertex : pi0",100,-2,2,100,-2,2);
	TH2F *hvpos4 = new TH2F("hvpos4","(x,y) projection of fitted decay vertex 2:Ds+",100,-2,2,100,-2,2);
	TH2F *hvpos5 = new TH2F("hvpos5","(x,y) projection of fitted decay vertex : Ds1*+",100,-2,2,100,-2,2);



	//
	// Now the analysis stuff comes...
	//
	
	
	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis();
	if (nevts==0) nevts= theAnalysis->GetEntries();
	
	// *** RhoCandLists for the analysis
	RhoCandList kplus, kminus, piplus, piminus, Dstar, Dsplist, Dsmlist, Ds2460, gammas, pi0list;
	
	// *** Mass selector for the jpsi cands

	RhoMassParticleSelector *DsMassSel=new RhoMassParticleSelector("D_s-",1.9686,0.5);
	RhoMassParticleSelector *DstarMassSel=new RhoMassParticleSelector("D_s*+", 2.1124,0.5);
	RhoMassParticleSelector *Ds2460MassSel=new RhoMassParticleSelector("D_s1+",2.457,0.3);
	RhoMassParticleSelector *DsMassSelPrefit=new RhoMassParticleSelector("D_s-",1.9686,0.5);
	RhoMassParticleSelector *DsMassSelPid=new RhoMassParticleSelector("D_s-",1.9686,0.3);
	RhoMassParticleSelector *DspMassSelPid=new RhoMassParticleSelector("D_s+",1.9686,0.3);
	RhoMassParticleSelector *Pi0MassSel=new RhoMassParticleSelector("pi0",0.134976,0.03);
	RhoEnergyParticleSelector *GamEnSel=new RhoEnergyParticleSelector("gam",0.05+5.,10.);
	RhoMomentumParticleSelector *Pi0MomSel=new RhoMomentumParticleSelector("pi0",0.1+5.,10.);

	// *** the lorentz vector of the initial Ds1+ + Ds- system
	TLorentzVector ini(0., 0.,9.58784,10.5719);
	Double_t kppx, kppy, kppz, kpE;
	Double_t kmpx, kmpy, kmpz, kmE;
	Double_t pmpx, pmpy, pmpz, pmE;
	Double_t m12, m23, m13, m2P12, m2P13, m2P23;
	
	// ***
	// the event loop
	// ***
	while (theAnalysis->GetEvent() && i++<nevts)
	{
		if ((i%100)==0) cout<<"evt " << i << endl;
				
		// *** Select with no PID info ('All'); type and mass are set 		
		theAnalysis->FillList(kplus,  "KaonAllPlus");
		theAnalysis->FillList(kminus, "KaonAllMinus");
		theAnalysis->FillList(piplus,  "PionAllPlus");
		theAnalysis->FillList(piminus, "PionAllMinus");
		theAnalysis->FillList(gammas, "Neutral");
		//cout <<"gam raw:"<<gammas.GetLength()<<endl;
		for (j=0;j<gammas.GetLength();++j) gam_en->Fill(gammas[j]->E());
		gammas.Select(GamEnSel);

		//cout <<"gam:"<<gammas.GetLength()<<endl;

		Dsplist.Cleanup();
		Dsmlist.Cleanup();
		// *** combinatorics for Ds-/+ to K+ K- pi-
		Dsplist.Combine(kminus,kplus,piplus);
	       	Dsmlist.Combine(kplus,kminus,piminus);

		pi0list.Cleanup();
		// *** combinatorics for pi0 to gg
		pi0list.Combine(gammas,gammas);
		pi0list.SetType(111);
		for (j=0;j<pi0list.GetLength();++j) pi0_mom->Fill(pi0list[j]->P());
		pi0list.Select(Pi0MomSel);
		//pi0list.Select(Pi0MassSel);
		
		Dstar.Cleanup();
		// *** combinatorics for Ds*+ to Ds+ gamma
		Dstar.Combine(Dsplist,gammas);
		
		Ds2460.Cleanup();
		// *** combinatorics for Ds1+ to Ds* K+
		Ds2460.Combine(Dstar,pi0list);

		// ***
		// *** do the TRUTH MATCH for pi0
		// ***
		gammas.SetType(22);
		gammas.Select(GamEnSel);
		//cout <<"pi0:"<<pi0list.GetLength()<<endl;
				
		for (j=0;j<pi0list.GetLength();++j) 
		{
			pi0_all->Fill( pi0list[j]->M() );

			if (theAnalysis->McTruthMatch(pi0list[j]))
			{ 
				pi0_ftm->Fill( pi0list[j]->M() );
				pi0_px_ftm->Fill(pi0list[j]->Px());
				pi0_py_ftm->Fill( pi0list[j]->Py());
				pi0_pz_ftm->Fill( pi0list[j]->Pz());
				pi0_E_ftm->Fill( pi0list[j]->E());
			 	pi0_mass_diff_ftm->Fill( pi0list[j]->GetMcTruth()->M() - pi0list[j]->M() );
				pi0_px_diff_ftm->Fill(pi0list[j]->GetMcTruth()->Px()- pi0list[j]->Px());
				pi0_py_diff_ftm->Fill( pi0list[j]->GetMcTruth()->Py()-pi0list[j]->Py());
				pi0_pz_diff_ftm->Fill( pi0list[j]->GetMcTruth()->Pz()-pi0list[j]->Pz());
				pi0_E_diff_ftm->Fill( pi0list[j]->GetMcTruth()->E() - pi0list[j]->E());
				pi0_Mom_ftm->Fill( pi0list[j]->P());
			}
			//			else 
			//	pi0_nm->Fill( pi0list[j]->M() );
			//	}
		}
		
		// *** do the TRUTH MATCH for Ds-
		// ***
		
		Dsmlist.SetType(-431);
		//cout <<"Ds-:"<<Dsmlist.GetLength()<<endl;
				
		for (j=0;j<Dsmlist.GetLength();++j) 
		{
			Dsm_all->Fill( Dsmlist[j]->M() );
			
			if (theAnalysis->McTruthMatch(Dsmlist[j]))
			{ 
				Dsm_ftm->Fill( Dsmlist[j]->M() );
			 	Dsm_mass_diff_ftm->Fill( Dsmlist[j]->GetMcTruth()->M() - Dsmlist[j]->M() );
				Dsm_px_ftm->Fill( Dsmlist[j]->Px());
				Dsm_py_ftm->Fill( Dsmlist[j]->Py());
				Dsm_pz_ftm->Fill( Dsmlist[j]->Pz());
				Dsm_E_ftm->Fill( Dsmlist[j]->E());
				Dsm_Mom_ftm->Fill( Dsmlist[j]->P());
				Dsm_px_diff_ftm->Fill(Dsmlist[j]->GetMcTruth()->Px() -  Dsmlist[j]->Px());
				Dsm_py_diff_ftm->Fill(Dsmlist[j]->GetMcTruth()->Py() -  Dsmlist[j]->Py());
				Dsm_pz_diff_ftm->Fill(Dsmlist[j]->GetMcTruth()->Pz() -  Dsmlist[j]->Pz());
				Dsm_E_diff_ftm->Fill(Dsmlist[j]->GetMcTruth()->E() -  Dsmlist[j]->E());

				RhoCandidate *kpluscand = Dsmlist[j]->Daughter(0);
				RhoCandidate *kminuscand = Dsmlist[j]->Daughter(1);
				RhoCandidate *piminuscand = Dsmlist[j]->Daughter(2);

				TLorentzVector m1(kppx, kppy, kppz, kpE); 
				TLorentzVector m2(kmpx, kmpy, kmpz, kmE); 
				TLorentzVector m3(pmpx, pmpy, pmpz, pmE);
				kppx = kpluscand->Px();
				kppy = kpluscand->Py();
				kppz = kpluscand->Pz();
				kpE = kpluscand->E();
				kmpx = kminuscand->Px();
				kmpy = kminuscand->Py();
				kmpz = kminuscand->Pz();
				kmE = kminuscand->E();
				pmpx = piminuscand->Px();
				pmpy = piminuscand->Py();
				pmpz = piminuscand->Pz();
				pmE = piminuscand->E();

				m12 = (m1+m2).M();
				m13 = (m1+m3).M();
				m23 = (m2+m3).M();

				m2P12 = m12*m12;
				m2P13 = m13*m13;
				m2P23 = m23*m23;

				h12_m->Fill(m12);
				h13_m->Fill(m13);
				h23_m->Fill(m23);

				hDalitzm12m13_m->Fill(m2P12,m2P13);
				hDalitzm12m23_m->Fill(m2P12,m2P23);
				hDalitzm13m23_m->Fill(m2P13,m2P23);
				
				hpmpx_m->Fill(pmpx);
				hpmpy_m->Fill(pmpy);
				hpmpz_m->Fill(pmpz);
				hpmE_m->Fill(pmE);
				
				hkmpx_m->Fill(kmpx);
				hkmpy_m->Fill(kmpy);
				hkmpz_m->Fill(kmpz);
				hkmE_m->Fill(kmE);
				
				hkppx_m->Fill(kppx);
				hkppy_m->Fill(kppy);
				hkppz_m->Fill(kppz);
				hkpE_m->Fill(kpE);
			
			}
			//	else 
			//			Dsm_nm->Fill( Dsmlist[j]->M() );
		}

		// *** do the TRUTH MATCH for Ds+
		// ***
		
		Dsplist.SetType(431);
		//cout <<"Dsp:"<<Dsplist.GetLength()<<endl;
				
		for (j=0;j<Dsplist.GetLength();++j) 
		{
			Dsp_all->Fill( Dsplist[j]->M() );
			
			if (theAnalysis->McTruthMatch(Dsplist[j]))
			{ 
				Dsp_ftm->Fill( Dsplist[j]->M() );
			 	Dsp_mass_diff_ftm->Fill( Dsplist[j]->GetMcTruth()->M() - Dsplist[j]->M() );
				Dsp_px_ftm->Fill( Dsplist[j]->Px());
				Dsp_py_ftm->Fill( Dsplist[j]->Py());
				Dsp_pz_ftm->Fill( Dsplist[j]->Pz());
				Dsp_E_ftm->Fill( Dsplist[j]->E());
				Dsp_Mom_ftm->Fill( Dsplist[j]->P());
				Dsp_px_diff_ftm->Fill(Dsplist[j]->GetMcTruth()->Px() -  Dsplist[j]->Px());
				Dsp_py_diff_ftm->Fill(Dsplist[j]->GetMcTruth()->Py() -  Dsplist[j]->Py());
				Dsp_pz_diff_ftm->Fill(Dsplist[j]->GetMcTruth()->Pz() -  Dsplist[j]->Pz());
				Dsp_E_diff_ftm->Fill(Dsplist[j]->GetMcTruth()->E() -  Dsplist[j]->E());

				RhoCandidate *kpluscand = Dsplist[j]->Daughter(0);
				RhoCandidate *kminuscand = Dsplist[j]->Daughter(1);
				RhoCandidate *pipluscand = Dsplist[j]->Daughter(2);

				TLorentzVector m1(kppx, kppy, kppz, kpE); 
				TLorentzVector m2(kmpx, kmpy, kmpz, kmE); 
				TLorentzVector m3(pmpx, pmpy, pmpz, pmE);
				kppx = kpluscand->Px();
				kppy = kpluscand->Py();
				kppz = kpluscand->Pz();
				kpE = kpluscand->E();
				kmpx = kminuscand->Px();
				kmpy = kminuscand->Py();
				kmpz = kminuscand->Pz();
				kmE = kminuscand->E();
				pmpx = pipluscand->Px();
				pmpy = pipluscand->Py();
				pmpz = pipluscand->Pz();
				pmE = pipluscand->E();

				m12 = (m1+m2).M();
				m13 = (m1+m3).M();
				m23 = (m2+m3).M();

				m2P12 = m12*m12;
				m2P13 = m13*m13;
				m2P23 = m23*m23;

				h12_p->Fill(m12);
				h13_p->Fill(m13);
				h23_p->Fill(m23);

				hDalitzm12m13_p->Fill(m2P12,m2P13);
				hDalitzm12m23_p->Fill(m2P12,m2P23);
				hDalitzm13m23_p->Fill(m2P13,m2P23);
				
				hpmpx_p->Fill(pmpx);
				hpmpy_p->Fill(pmpy);
				hpmpz_p->Fill(pmpz);
				hpmE_p->Fill(pmE);
				
				hkmpx_p->Fill(kmpx);
				hkmpy_p->Fill(kmpy);
				hkmpz_p->Fill(kmpz);
				hkmE_p->Fill(kmE);
				
				hkppx_p->Fill(kppx);
				hkppy_p->Fill(kppy);
				hkppz_p->Fill(kppz);
				hkpE_p->Fill(kpE);
			
			}
			//		else 
			//	Dsp_nm->Fill( Dsplist[j]->M() );
		}



			
		// ***
		// *** do the TRUTH MATCH for Ds*+
		// ***

		
		Dstar.SetType(433);
		pi0list.SetType(111);
		gammas.SetType(22);
				
		//cout <<"D*:"<<Dstar.GetLength()<<endl;
		
		for (j=0;j<Dstar.GetLength();++j) 
		{
			Dstar_all->Fill( Dstar[j]->M() );

			if (theAnalysis->McTruthMatch(Dstar[j]))
			{ 
				Dstar_ftm->Fill( Dstar[j]->M() );
			 	Dstar_mass_diff_ftm->Fill( Dstar[j]->GetMcTruth()->M() - Dstar[j]->M() );
				Dstar_px_ftm->Fill( Dstar[j]->Px());
				Dstar_py_ftm->Fill( Dstar[j]->Py());
				Dstar_pz_ftm->Fill( Dstar[j]->Pz());
				Dstar_E_ftm->Fill( Dstar[j]->E());
				Dstar_Mom_ftm->Fill( Dstar[j]->P());
				Dstar_px_diff_ftm->Fill( Dstar[j]->GetMcTruth()->Px() - Dstar[j]->Px());
				Dstar_py_diff_ftm->Fill( Dstar[j]->GetMcTruth()->Py() - Dstar[j]->Py());
				Dstar_pz_diff_ftm->Fill( Dstar[j]->GetMcTruth()->Pz() - Dstar[j]->Pz());
				Dstar_E_diff_ftm->Fill( Dstar[j]->GetMcTruth()->E() - Dstar[j]->E());
				
			}
			//	else 
			//	Dstar_nm->Fill( Dstar[j]->M() );
		}
	// ***

		
		// ***
		// *** do the TRUTH MATCH for Ds'1+
		// ***
		
		Ds2460.SetType(10433);
		//cout<<Ds2460.PdgCode()<<endl;
		//cout <<"Ds2460:"<<Ds2460.GetLength()<<endl;
						
		for (j=0;j<Ds2460.GetLength();++j) 
		{
			D2460_all->Fill( Ds2460[j]->M() );
			
			if (theAnalysis->McTruthMatch(Ds2460[j]))
			{ 
				D2460_ftm->Fill( Ds2460[j]->M() );
				D2460_mass_diff_ftm->Fill( Ds2460[j]->GetMcTruth()->M() - Ds2460[j]->M() );
				D2460_px_ftm->Fill( Ds2460[j]->Px());
				D2460_py_ftm->Fill( Ds2460[j]->Py());
				D2460_pz_ftm->Fill( Ds2460[j]->Pz());
				D2460_E_ftm->Fill( Ds2460[j]->E());
				D2460_Mom_ftm->Fill( Ds2460[j]->P());
				D2460_px_diff_ftm->Fill( Ds2460[j]->GetMcTruth()->Px() - Ds2460[j]->Px());
				D2460_py_diff_ftm->Fill( Ds2460[j]->GetMcTruth()->Py() - Ds2460[j]->Py());
				D2460_pz_diff_ftm->Fill( Ds2460[j]->GetMcTruth()->Pz() - Ds2460[j]->Pz());
				D2460_E_diff_ftm->Fill( Ds2460[j]->GetMcTruth()->E() - Ds2460[j]->E());	
				//	hvpos->Fill(Ds2460[j]->Pos()->X(),Ds2460[j]->Pos()->Y());
			
			}
			//	else 
			//	D2460_nm->Fill( Ds2460[j]->M() );
		}
		

		Dsmlist.Select(DsMassSelPid);
		Dsplist.Select(DsMassSelPid);
		Dstar.Select(DstarMassSel);
		pi0list.Select(Pi0MassSel);
		gammas.Select(GamEnSel);
		pi0list.Select(Pi0MomSel);

		
		for (j=0;j<Dsmlist.GetLength();++j) 
		{
			PndKinVtxFitter vtxfitter(Dsmlist[j]);	// instantiate a vertex fitter
			vtxfitter.Fit();
			
			double chi2_vtx=vtxfitter.GetChi2();	// access chi2 of fit
			double prob_vtx=vtxfitter.GetProb();
			Dsm_chi2_vf->Fill(chi2_vtx);
			Dsm_chi2_prob_vf->Fill(prob_vtx);
		
			
			if (chi2_vtx>0.&&chi2_vtx<1.)				// when good enough, fill some histos
			{
				RhoCandidate *jfit = Dsmlist[j]->GetFit();	// access the fitted cand
				
				TVector3 jVtx=jfit->Pos();		// and the decay vertex position
				
				Dsm_vf->Fill(jfit->M());    
				
				hvpos2m->Fill(jVtx.X(),jVtx.Y());
				
				Dsm_vtx_z->Fill(jVtx.Z());
				
				RhoCandidate *mcdsm = Dsmlist[j]->GetMcTruth();
				if (mcdsm && mcdsm->PdgCode()==-431)
				  {hm_diff2m->Fill( mcdsm->M() - jfit->M());
				  }
				
			}
			
			
		}
		/*
		// ***
		// *** do VERTEX FIT (Ds1+)

		// ***

		//Ds2460.Cleanup();
		for (j=0;j<Ds2460.GetLength();++j) 
		{
			PndKinVtxFitter vtxfitter(Ds2460[j]);	// instantiate a vertex fitter
			vtxfitter.Fit();
			
			double chi2_vtx=vtxfitter.GetChi2();	// access chi2 of fit
			double prob_vtx=vtxfitter.GetProb();
			D2460_chi2_vf->Fill(chi2_vtx);
			D2460_chi2_prob_vf->Fill(prob_vtx);
			if (chi2_vtx>0.&&chi2_vtx<1.)				// when good enough, fill some histos
			{
				RhoCandidate *jfit = Ds2460[j]->GetFit();	// access the fitted cand
				TVector3 jVtx=jfit->Pos();		// and the decay vertex position
				
				D2460_vf->Fill(jfit->M());            
				hvpos5->Fill(jVtx.X(),jVtx.Y());
				D2460_vf->Fill(jfit->Z());
			RhoCandidate *mcd = D2460[j]->GetMcTruth();
			if (mcd && mcd->PdgCode()==10433){
			  hm_diff3->Fill( mcd->M() - jVtx->M());
			  }

			}
		}

		*/


		// ***
		// *** do MASS CONSTRAINT FIT (Ds-/+)

		
		// ***
		//Dsmlist.Cleanup();
		for (j=0;j<Dsmlist.GetLength();++j) 
		{
			PndKinFitter mfitter(Dsmlist[j]);        // instantiate the PndKinFitter in pp bar
			mfitter.AddMassConstraint(1.9686);
			// add the mass constraint
			mfitter.Fit();				// do fit
		
			double chi2_m = mfitter.GetChi2();	// get chi2 of fit
			Dsm_chi2_mf->Fill(chi2_m);
			double prob_m = mfitter.GetProb();
						
			if (chi2_m>0.&&chi2_m<1)				// when good enough, fill some histo
			{
				RhoCandidate *jfit = Dsmlist[j]->GetFit();// access the fitted cand
				Dsm_chi2_prob_mf->Fill(prob_m);
				if(prob_m>0.01){Dsm_mcf->Fill(jfit->M());}
			}

			RhoCandidate *mcdsm = Dsmlist[j]->GetMcTruth();
			if (mcdsm && mcdsm->PdgCode()==-431)
			  {hm_diff2m_mcf->Fill( mcdsm->M() - jfit->M());
			  }
		}
		 
			for (j=0;j<Dsplist.GetLength();++j) 
		{
			PndKinFitter mfitter(Dsplist[j]);        // instantiate the PndKinFitter in pp bar
			mfitter.AddMassConstraint(1.9686);
			// add the mass constraint
			mfitter.Fit();				// do fit
		
			double chi2_m = mfitter.GetChi2();	// get chi2 of fit
			Dsp_chi2_mf->Fill(chi2_m);
			double prob_m = mfitter.GetProb();
						
			if (chi2_m>0.&&chi2_m<1)				// when good enough, fill some histo
			{
				RhoCandidate *jfit = Dsplist[j]->GetFit();// access the fitted cand
				Dsp_chi2_prob_mf->Fill(prob_m);
				if(prob_m>0.01){Dsp_mcf->Fill(jfit->M());}
			}
			RhoCandidate *mcdsm = Dsplist[j]->GetMcTruth();
			if (mcdsm && mcdsm->PdgCode()==431)
			  {hm_diff2p_mcf->Fill( mcdsm->M() - jfit->M());
			  }

		}
		   
		

	       
		// *** do MASS CONSTRAINT FIT (Ds*+)
		// ***
		//	Dstar.Cleanup();
		for (j=0;j<Dstar.GetLength();++j) 
		{
			PndKinFitter mfitter(Dstar[j]);        // instantiate the PndKinFitter in pp bar
			mfitter.AddMassConstraint(2.1124);
			// add the mass constraint
			mfitter.Fit();				// do fit
		
			double chi2_m = mfitter.GetChi2();	// get chi2 of fit
			Dstar_chi2_mf->Fill(chi2_m);
			double prob_m = mfitter.GetProb();
			if (chi2_m>0.&&chi2_m<1.)				// when good enough, fill some histo
			{
				RhoCandidate *jfit = Dstar[j]->GetFit();// access the fitted cand
				Dstar_chi2_prob_mf->Fill(prob_m);
		
				TVector3 jVtx=jfit->Pos();
				if(prob_m>0.01){Dstar_mcf->Fill(jfit->M());}

			}
			RhoCandidate *mcdsm = Dstar[j]->GetMcTruth();
			if (mcdsm && mcdsm->PdgCode()==433)
			  {hm_diffstar_mcf->Fill( mcdsm->M() - jfit->M());
			  }
		}
		
	       
		// *** do MASS CONSTRAINT FIT (pi0)
		// ***
		//	Dstar.Cleanup();
		for (j=0;j<pi0list.GetLength();++j) 
		{
			PndKinFitter mfitter(pi0list[j]);        // instantiate the PndKinFitter in pp bar
			mfitter.AddMassConstraint(0.134976);
			// add the mass constraint
			mfitter.Fit();				// do fit
		
			double chi2_m = mfitter.GetChi2();	// get chi2 of fit
			pi0_chi2_mf->Fill(chi2_m);
			double prob_m = mfitter.GetProb();
			if (chi2_m>0.&&chi2_m<1.)				// when good enough, fill some histo
			{
				RhoCandidate *jfit = pi0list[j]->GetFit();// access the fitted cand
				pi0_chi2_prob_mf->Fill(prob_m);
		
				TVector3 jVtx=jfit->Pos();
				if(prob_m>0.01){pi0_mcf->Fill(jfit->M());}

			}

			RhoCandidate *mcpi0 = pi0list[j]->GetMcTruth();
			if (mcpi0 && mcpi0->PdgCode()==111)
			  {hm_diffpi0_mcf->Fill( mcpi0->M() - jfit->M());
			  }
		}
		

		//Ds2460.Cleanup();
		Ds2460.Select(Ds2460MassSel);
		// *** do MASS CONSTRAINT FIT (D'1s+)
		// ***
		for (j=0;j<Ds2460.GetLength();++j) 
		{
			PndKinFitter mfitter(Ds2460[j]);        // instantiate the PndKinFitter in pp bar
			mfitter.AddMassConstraint(2.457);
			// add the mass constraint
			mfitter.Fit();				// do fit
		
			double chi2_m = mfitter.GetChi2();	// get chi2 of fit
			D2460_chi2_mf->Fill(chi2_m);
			double prob_m = mfitter.GetProb();
			if (chi2_m>0&&chi2_m<1)				// when good enough, fill some histo
			{
				RhoCandidate *jfit = Ds2460[j]->GetFit();// access the fitted cand
				
				D2460_chi2_prob_mf->Fill(prob_m);
				if(prob_m>0.01){D2460_mcf->Fill(jfit->M());}
			
			}
			RhoCandidate *mcdsm = Ds2460[j]->GetMcTruth();
			if (mcdsm && mcdsm->PdgCode()==10433)
			  {hm_diffds1_mcf->Fill( mcdsm->M() - jfit->M());
			  }
		}




		// ***
		// *** TRUE PID combinatorics
		// ***
				
		// *** do MC truth match for PID type
		SelectTruePid(theAnalysis, kplus);
		SelectTruePid(theAnalysis, kminus);
		SelectTruePid(theAnalysis, piplus);
		SelectTruePid(theAnalysis, piminus);
		SelectTruePid(theAnalysis, gammas);
			

		gammas.Select(GamEnSel);
		pi0list.Select(Pi0MomSel);
		pi0list.Select(Pi0MassSel);
		// *** all combinatorics again with true PID
		pi0list.Combine(gammas, gammas);
		for (j=0;j<pi0list.GetLength();++j){ 
		  if(gammas.GetLength()>50) continue;
		    pi0_trpid->Fill( pi0list[j]->M() );
		}
		Dsmlist.Combine(kplus, kminus, piminus);
		for (j=0;j<Dsmlist.GetLength();++j) Dsm_trpid->Fill( Dsmlist[j]->M() );
		Dsmlist.Select(DsMassSelPid);

		Dsplist.Combine(kplus, kminus, piplus);
		for (j=0;j<Dsplist.GetLength();++j) Dsp_trpid->Fill( Dsplist[j]->M() );
		Dsplist.Select(DsMassSelPid);
	
		Dstar.Combine(Dsplist,gammas);
		for (j=0;j<Dstar.GetLength();++j){ 
		  if(gammas.GetLength()>50) continue;
		    Dstar_trpid->Fill( Dstar[j]->M() );
		    Dstar.Select(DstarMassSel);
		}
		Ds2460.Combine(Dstar, pi0list);
		for (j=0;j<Ds2460.GetLength();++j){ 
		  if(pi0list.GetLength()>50) continue;
		    D2460_trpid->Fill( Ds2460[j]->M() );
		  
		}




		// ***
		// *** BEST PID combinatorics
		// ***
		
		// *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection
		theAnalysis->FillList(kplus,  "KaonBestPlus",  "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(kminus, "KaonBestMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(piplus,  "PionBestPlus",  "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(piminus, "PionBestMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(gammas, "Neutral", "PidAlgoEmcBayes");

		pi0list.Combine(gammas,gammas);
		pi0list.Select(Pi0MassSel);
		pi0list.Select(Pi0MomSel);
		Dsmlist.Select(DspMassSelPid);
		Dsplist.Select(DsMassSelPid);
		Dstar.Select(DstarMassSel);
		for (j=0;j<pi0list.GetLength();++j){
		  if(gammas.GetLength()>50) continue;
		  pi0_best->Fill( pi0list[j]->M() );
		}
		
		Dsmlist.Combine(kplus, kminus, piminus);
		for (j=0;j<Dsmlist.GetLength();++j){
		  // Dsmlist.Select(DspMassSelPid);
		  Dsm_best->Fill( Dsmlist[j]->M() );
		  
		}
		Dsplist.Combine(kplus, kminus, piplus);
		for (j=0;j<Dsplist.GetLength();++j){
		  //Dsplist.Select(DsMassSelPid);
		  Dsp_best->Fill( Dsplist[j]->M() );
		  
		}
		Dstar.Combine(Dsplist, gammas);
		for (j=0;j<Dstar.GetLength();++j){
		  if(gammas.GetLength()>50) continue;
		  //Dstar.Select(DstarMassSel);
		  Dstar_best->Fill( Dstar[j]->M() );
		  
		}
	
		Ds2460.Combine(Dstar, pi0list);
		for (j=0;j<Ds2460.GetLength();++j){
		  if(gammas.GetLength() > 50) continue;
		  if(pi0list.GetLength() > 50) continue;
		  D2460_best->Fill( Ds2460[j]->M() );
		}
	}
	


	// *** write out all the histos

	out->cd();
	

	D2460_all->Write();
	Dsm_all->Write();
	Dsp_all->Write();
	Dstar_all->Write();
	pi0_all->Write();
	
	D2460_best->Write();
	Dsm_best->Write();
	Dsp_best->Write();
	pi0_best->Write();
	Dstar_best->Write();

	D2460_trpid->Write();
	Dsm_trpid->Write();
	Dsp_trpid->Write();
	pi0_trpid->Write();
	Dstar_trpid->Write();

	D2460_ftm->Write();
	Dsp_ftm->Write();
	Dsm_ftm->Write();
	pi0_ftm->Write();
	Dstar_ftm->Write();
	//	hm_diff->Write();
	hm_diff2m->Write();
	hm_diff2m_mcf->Write();
	hm_diff2p_mcf->Write();
	hm_diffpi0_mcf->Write();
	hm_diffstar_mcf->Write();
	hm_diffds1_mcf->Write();
	//	hm_diff2p->Write();
	//	hm_diff3->Write();
	//	hm_diff4->Write();


	D2460_mcf->Write();
	//	D2460_4cf->Write();
	D2460_chi2_mf->Write();
	//	D2460_chi2_4c->Write();
	//	D2460_chi2_vf->Write();
	//D2460_chi2_prob_4c->Write();
	D2460_chi2_prob_mf->Write();
	
	//	pi0_vf->Write();
	pi0_mcf->Write();
	pi0_chi2_mf->Write();
	//	pi0_chi2_vf->Write();
	//	pi0_chi2_prob_vf->Write();
	pi0_chi2_prob_mf->Write();
	
	Dstar_mcf->Write();
	Dstar_chi2_mf->Write();
	//Dstar_chi2_vf->Write();
	Dstar_chi2_prob_mf->Write();

	Dsm_mcf->Write();
	//	Dsm_vf->Write();
	Dsm_chi2_mf->Write();
	//	Dsm_chi2_vf->Write();
	Dsm_chi2_prob_mf->Write();
	//Dsm_chi2_prob_vf->Write();

	Dsp_mcf->Write();
	//Dsp_vf->Write();
	Dsp_chi2_mf->Write();
	//Dsp_chi2_vf->Write();
	Dsp_chi2_prob_mf->Write();
	//Dsp_chi2_prob_vf->Write();



	hvpos2m->Write();


	//	pi0_vtx_z->Write();
	//Dsp_vtx_z->Write();
	Dsm_vtx_z->Write();
	//Dstar_vtx_z->Write();
	//D2460_vtx_z->Write();
	
	D2460_px_ftm->Write();
	D2460_py_ftm->Write();
	D2460_pz_ftm->Write();
	D2460_E_ftm->Write();
	D2460_px_diff_ftm->Write();
	D2460_py_diff_ftm->Write();
	D2460_pz_diff_ftm->Write();
	D2460_E_diff_ftm->Write();
	D2460_mass_diff_ftm->Write();

	pi0_px_ftm->Write();
	pi0_py_ftm->Write();
	pi0_pz_ftm->Write();
	pi0_E_ftm->Write();
	pi0_px_diff_ftm->Write();
	pi0_py_diff_ftm->Write();
	pi0_pz_diff_ftm->Write();
	pi0_E_diff_ftm->Write();
	pi0_mass_diff_ftm->Write();
	
	
	Dstar_px_ftm->Write();
	Dstar_py_ftm->Write();
	Dstar_pz_ftm->Write();
	Dstar_E_ftm->Write();
	Dstar_px_diff_ftm->Write();
	Dstar_py_diff_ftm->Write();
	Dstar_pz_diff_ftm->Write();
	Dstar_E_diff_ftm->Write();
	Dstar_mass_diff_ftm->Write();
	
	Dsm_px_ftm->Write();
	Dsm_py_ftm->Write();
	Dsm_pz_ftm->Write();
	Dsm_E_ftm->Write();
	Dsm_px_diff_ftm->Write();
	Dsm_py_diff_ftm->Write();
	Dsm_pz_diff_ftm->Write();
	Dsm_E_diff_ftm->Write();
	Dsm_mass_diff_ftm->Write();
	Dsp_px_ftm->Write();
	Dsp_py_ftm->Write();
	Dsp_pz_ftm->Write();
	Dsp_E_ftm->Write();
	Dsp_px_diff_ftm->Write();
	Dsp_py_diff_ftm->Write();
	Dsp_pz_diff_ftm->Write();
	Dsp_E_diff_ftm->Write();
	Dsp_mass_diff_ftm->Write();


	hkmpx_m->Write();
	hkmpy_m->Write();
	hkmpz_m->Write();
	hkmE_m->Write();
	hkppx_m->Write();
	hkppy_m->Write();
	hkppz_m->Write();
	hkpE_m->Write();
	hpmpx_m->Write();
	hpmpy_m->Write();
	hpmpz_m->Write();
	hpmE_m->Write();

	hkmpx_p->Write();
	hkmpy_p->Write();
	hkmpz_p->Write();
	hkmE_p->Write();
	hkppx_p->Write();
	hkppy_p->Write();
	hkppz_p->Write();
	hkpE_p->Write();
	hpmpx_p->Write();
	hpmpy_p->Write();
	hpmpz_p->Write();
	hpmE_p->Write();

	h12_m->Write();
	h13_m->Write();
	h23_m->Write();

	h12_p->Write();
	h13_p->Write();
	h23_p->Write();

	hDalitzm12m13_m->Write();
	hDalitzm12m23_m->Write();
	hDalitzm13m23_m->Write();

	hDalitzm12m13_p->Write();
	hDalitzm12m23_p->Write();
	hDalitzm13m23_p->Write();


	pi0_mom->Write();
	gam_en->Write();


	out->Save();

}