////////////////////////////////////////////////////////////////////////
//                                                                    // 
//Elisabetta 15/08/2013                                               //
//macro modified adding a root tree, respect to the standard tutorial //
//Analysis: p pbar to Ds1+ Ds-                                        //
//where: Ds1(2317) ---> 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_anaDs2317(int nevts)
{
	// *** 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_anaDs2317.root","RECREATE");
		
	// *** create some histograms



	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 *D2317_lpid = new TH1F("D2317_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 *D2317_best = new TH1F("D2317_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);

	//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,-1.5,1.5);
	TH1F *pi0_E_ftm  = new TH1F("pi0_E_ftm","#pi0 E ",100,0,10);
	

	//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);


	
	//Ds2317 plots: all and ftm
	TH1F *D2317_all  = new TH1F("D2317_all","Ds'1^{+} mass (all)",100,0.,5.);
	TH1F *D2317_px_ftm  = new TH1F("D2317_px_ftm","Ds'1+ px ",100,-1.5,1.5);
	TH1F *D2317_py_ftm  = new TH1F("D2317_py_ftm","Ds'1+ py ",100,-1.5,1.5);
	TH1F *D2317_pz_ftm  = new TH1F("D2317_pz_ftm","Ds'1+ pz ",100,-1,9);
	TH1F *D2317_E_ftm  = new TH1F("D2317_E_ftm","Ds'1+ E ",100,0,10);
	TH1F *D2317_mass_diff_ftm  = new TH1F("D2317_mass_diff_ftm","Ds'1+ mass resolution (full truth match)",100,-0.2,0.2);
	TH1F *D2317_ftm  = new TH1F("D2317_ftm","Ds'1+ mass (full truth match)",100,0.,4.);
	TH1F *D2317_px_diff_ftm  = new TH1F("D2317_px_diff_ftm","Ds'1+ px resolution",100,-0.5,0.5);
	TH1F *D2317_py_diff_ftm  = new TH1F("D2317_py_diff_ftm","Ds'1+ py resolution",100,-0.5,0.5);
	TH1F *D2317_pz_diff_ftm  = new TH1F("D2317_pz_diff_ftm","Ds'1+ pz resolution",100,-0.5,0.5);
	TH1F *D2317_E_diff_ftm  = new TH1F("D2317_E_diff_ftm","Ds'1+ E resolution",100,-0.5,0.5);
	TH1F *D2317_diff_mcf  = new TH1F("D2317_diff_mcf","Ds'1+ mass resolution",100,-0.5,0.5);

	//best 

	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 *D2317_tpid  = new TH1F("D2317_tpid"," Ds1^{+} mass (tight pid)",100,0,5);
	
	TH1F *D2317_trpid = new TH1F("D2317_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 *D2317_nm  = new TH1F("D2317_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);

	//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 *D2317_mfc   = new TH1F("D2317_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 *D2317_4cf  = new TH1F("D2317_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_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,3);
	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 *D2317_mcf  = new TH1F("D2317_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 *D2317_chi2_mf  = new TH1F("D2317_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 *D2317_chi2_vf  = new TH1F("D2317_chi2_vf", "D2317: #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 *D2317_vtx_z   = new TH1F("D2317_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 *D2317_chi2_prob_vf   = new TH1F("D2317_chi2_prob_vf",  "D's1+: prob vertex fit",100,-0.01,1);
	TH1F *D2317_chi2_prob_mf   = new TH1F("D2317_chi2_prob_mf",  "Ds'1+: prob mass fit",100,-0.01,1);
	TH1F *D2317_chi2_4c   = new TH1F("D2317_chi2_4c",  "Ds'1+: #chi^{2} 4C fit",100,-1,20);
	TH1F *D2317_chi2_prob_4c   = new TH1F("D2317_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, Ds2317, gammas, pi0list;
	
	// *** Mass selector for the jpsi cands

	RhoMassParticleSelector *DsMassSel=new RhoMassParticleSelector("D_s-",1.9686,0.5);
	RhoMassParticleSelector *Ds2317MassSel=new RhoMassParticleSelector("D_s0+",2.317,0.4);
	RhoMassParticleSelector *DsMassSelPrefit=new RhoMassParticleSelector("D_s-",1.9686,0.5);
	RhoMassParticleSelector *DsMassSelPid=new RhoMassParticleSelector("D_s-",1.9686,1.0);
	RhoMassParticleSelector *Pi0MassSel=new RhoMassParticleSelector("pi0",0.134976,0.1);

	

	//il pi0 e' largo <2MeV; la Ds molto meno.
	// *** 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;
	Double_t energy_g;
	Int_t counting_g;
	// ***
	// the event loop
	// ***
	while (theAnalysis->GetEvent() && i++<nevts)
	  {
	    if ((i%10)==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");
	    
	    
	    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);
	
	
	    Ds2317.Cleanup();
	    // *** combinatorics for Ds1+ to Ds* K+
	    Ds2317.Combine(Dsplist,pi0list);
	    
	    // ***
	    // *** do the TRUTH MATCH for pi0
	    // ***
	    //pi0list.Cleanup();
	    pi0list.SetType(111);
	    gammas.SetType(22);
	    
	    for (j=0;j<gammas.GetLength();++j) 
	      {
		energy_g = gammas[j]->E();
		counting_g = gammas.GetLength();
		if(energy_g<0.15) continue;
		if(counting_g>100) continue;
	      }
	
	    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());
		    
		    
		  }
		
	      }
	  
	    // *** do the TRUTH MATCH for Ds-
		// ***
		//Dsmlist.Cleanup();
	    Dsmlist.SetType(-431);
	    
	    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_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);
		    
		  }
	      }	  
		// *** do the TRUTH MATCH for Ds+
		// ***
		//Dsplist.Cleanup();
	    Dsplist.SetType(431);
	    
	    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_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);
		    
		  }

	      }
	    
	  


		
		// ***
		// *** do the TRUTH MATCH for Ds'1+
		// ***
		//Ds2317.Cleanup();
	    Ds2317.SetType(10431);
	    //cout<<Ds2317.PdgCode()<<endl;
	    
	    for (j=0;j<Ds2317.GetLength();++j) 
	      {
		D2317_all->Fill( Ds2317[j]->M() );
		
		if (theAnalysis->McTruthMatch(Ds2317[j]))
		  { 
		    D2317_ftm->Fill( Ds2317[j]->M() );
		    D2317_mass_diff_ftm->Fill( Ds2317[j]->GetMcTruth()->M() - Ds2317[j]->M() );
		    D2317_px_ftm->Fill( Ds2317[j]->Px());
		    D2317_py_ftm->Fill( Ds2317[j]->Py());
		    D2317_pz_ftm->Fill( Ds2317[j]->Pz());
		    D2317_E_ftm->Fill( Ds2317[j]->E());
		    D2317_px_diff_ftm->Fill( Ds2317[j]->GetMcTruth()->Px() - Ds2317[j]->Px());
		    D2317_py_diff_ftm->Fill( Ds2317[j]->GetMcTruth()->Py() - Ds2317[j]->Py());
		    D2317_pz_diff_ftm->Fill( Ds2317[j]->GetMcTruth()->Pz() - Ds2317[j]->Pz());
		    D2317_E_diff_ftm->Fill( Ds2317[j]->GetMcTruth()->E() - Ds2317[j]->E());	
		    //	hvpos->Fill(Ds2317[j]->Pos()->X(),Ds2317[j]->Pos()->Y());
		    
		  }


	      }
	    //
	    // for D*s an Ds- I do mass constraint
	    // D's+1 for now got only vertex fit
	    //
	    //
	    
	    Dsplist.Select(DsMassSelPrefit);
	    Dsmlist.Select(DsMassSelPrefit);
	    pi0list.Select(Pi0MassSel);
	    
	    

		// ***
		// *** do VERTEX FIT (Ds1+)

		// ***
	
	    
		//Ds2317.Cleanup();
	    for (j=0;j<Ds2317.GetLength();++j) 
	      {
		PndKinVtxFitter vtxfitter(Ds2317[j]);	// instantiate a vertex fitter
		vtxfitter.Fit();
		
		double chi2_vtx=vtxfitter.GetChi2();	// access chi2 of fit
		double prob_vtx=vtxfitter.GetProb();
		D2317_chi2_vf->Fill(chi2_vtx);
		D2317_chi2_prob_vf->Fill(prob_vtx);
		if (chi2_vtx>0.&&chi2_vtx<2.)				// when good enough, fill some histos
		  {
		    RhoCandidate *jfit = Ds2317[j]->GetFit();	// access the fitted cand
		    TVector3 jVtx=jfit->Pos();		// and the decay vertex position
		    
		    D2317_vf->Fill(jfit->M());            
		    hvpos5->Fill(jVtx.X(),jVtx.Y());
		    D2317_vf->Fill(jfit->Z());
		    RhoCandidate *mcd = D2317[j]->GetMcTruth();
		    if (mcd && mcd->PdgCode()==10431){
		      hm_diff3->Fill( mcd->M() - jVtx->M());
		    }
		  }
	      }
	      
	        
	    
	    // *** some rough mass selection
	    Dsmlist.Select(DsMassSel);
	    Dsplist.Select(DsMassSel);
	    pi0list.Select(Pi0MassSel);

			
		
		//Ds2317.Cleanup();
	    Ds2317.Select(Ds2317MassSel);
		// *** do MASS CONSTRAINT FIT (D'1s+)
		// ***
	    for (j=0;j<Ds2317.GetLength();++j) 
	      {
		PndKinFitter mfitter(Ds2317[j]);        // instantiate the PndKinFitter in pp bar
		mfitter.AddMassConstraint(2.317);
		// add the mass constraint
		mfitter.Fit();				// do fit
		
		double chi2_m = mfitter.GetChi2();	// get chi2 of fit
		D2317_chi2_mf->Fill(chi2_m);
		double prob_m = mfitter.GetProb();
		if (chi2_m>0&&chi2_m<2)				// when good enough, fill some histo
		  {
		    RhoCandidate *jfit = Ds2317[j]->GetFit();// access the fitted cand
		    
		    D2317_chi2_prob_mf->Fill(prob_m);
		    if(prob_m>0.01){D2317_mcf->Fill(jfit->M());}
		    
		  }
	      }
	    
	    
	    Dsmlist.Select(DsMassSelPid);
	    Dsplist.Select(DsMassSelPid);
	    pi0list.Select(Pi0MassSel);
	    
	    // ***
	    // *** 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);
	    
	    // *** all combinatorics again with true PID
	    pi0list.Combine(gammas, gammas);
	    for (j=0;j<pi0list.GetLength();++j){
	      if(gammas.GetLength()<50&&pi0list[j]->Daughter(0)->P()>0.150&&pi0list[j]->Daughter(1)->P()>0.150){
		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);
	    
	    Ds2317.Combine(Dsplist, pi0list);
	    if(pi0list.GetLength()<50){
	      for (j=0;j<Ds2317.GetLength();++j) D2317_trpid->Fill( Ds2317[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);
	      for (j=0;j<pi0list.GetLength();++j){
		//	int mydau = pi0list[j]->NDaughters();
		//int mylist = gammas.GetLength();
		//	if(mylist<50){
		  pi0_best->Fill( pi0list[j]->M() );
		  //}
	      }
	      Dsmlist.Combine(kplus, kminus, piminus);
	      for (j=0;j<Dsmlist.GetLength();++j){
		//int mydau = Dsmlist[j]->NDaughters();
		//int mylistdsm = Dsmlist.GetLength();
		//if(mylistdsm<200&&mydau==3){
		  
		  Dsm_best->Fill( Dsmlist[j]->M() );
		  Dsmlist.Select(DsMassSelPid);
		  //	}
	      }
	      Dsplist.Combine(kplus, kminus, piplus);
	      for (j=0;j<Dsplist.GetLength();++j){
		//int mydau = Dsplist[j]->NDaughters();
		///int mylistdsm = Dsplist.GetLength();
		//if(mylistdsm<200&&mydau==3){ 
		  
		  Dsp_best->Fill( Dsplist[j]->M() );
		  Dsplist.Select(DsMassSelPid);	
		}
	      // }
	      
	      Ds2317.Combine(Dsplist, pi0list);
	      for (j=0;j<Ds2317.GetLength();++j){
		//	int mydau = Ds2317[j]->NDaughters();
		//int mylistpi0 = pi0list.GetLength();
		//int dsp = Dsplist.GetLength();
		//if(mylistpi0<100&&mydau==2){
		  
		  D2317_best->Fill( Ds2317[j]->M() );
		  //	}
	      }
	      
	  
	  }
	    // *** write out all the histos
	out->cd();
	
	D2317_all->Write();
	Dsm_all->Write();
	Dsp_all->Write();
	pi0_all->Write();
	D2317_best->Write();
	Dsm_best->Write();
	Dsp_best->Write();
	pi0_best->Write();

	D2317_trpid->Write();
	Dsm_trpid->Write();
	Dsp_trpid->Write();
	pi0_trpid->Write();

	D2317_ftm->Write();
	Dsp_ftm->Write();
	Dsm_ftm->Write();
	pi0_ftm->Write();

	D2317_mcf->Write();
	D2317_chi2_mf->Write();
	D2317_chi2_prob_mf->Write();
	
	D2317_vtx_z->Write();
	
	D2317_px_ftm->Write();
	D2317_py_ftm->Write();
	D2317_pz_ftm->Write();
	D2317_E_ftm->Write();
	D2317_px_diff_ftm->Write();
	D2317_py_diff_ftm->Write();
	D2317_pz_diff_ftm->Write();
	D2317_E_diff_ftm->Write();
	D2317_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();
	

	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();


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

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


	out->Save();
	//	out->Close();
}