//class PndAnalysis;
//class PndAnaPidSelector;
//class RhoCandList;
//class RhoTuple;

#include "PndTutAnaTask.h"

// C++ headers
#include <string>
#include <iostream>
#include <map>

//Fair headers
#include "FairLogger.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairParRootFileIo.h"
#include "FairParAsciiFileIo.h"

//ROOT headers
#include "TFile.h"
#include "TLorentzVector.h"
#include "TDatabasePDG.h"
#include "TVector3.h"

//Rho headers
#include "RhoHistogram/RhoTuple.h"
#include "PndRhoTupleQA.h"
#include "RhoCandList.h"
#include "RhoCandidate.h"
#include "RhoMassParticleSelector.h"
#include "PndVtxPoca.h"

//Analysis headers
#include "PndAnalysis.h"
#include "PndEventShape.h"
#include "PndKinVtxFitter.h"
#include "PndKalmanVtxFitter.h"
#include "PndKinFitter.h"
#include "PndPidCandidate.h"
//#include "PndTrack.h"


using std::cout;
using std::endl;

ClassImp(PndTutAnaTask)

//---------------Default constructor------------
PndTutAnaTask::PndTutAnaTask():
	FairTask("Panda Tutorial Analysis Task") {
	}
//------------------------------------------

//--------------- destructor ----------------
PndTutAnaTask::~PndTutAnaTask(){}
//-----------------------------------------


//enum pidNumbers {
//
//
//
//
//
//};

void PndTutAnaTask::CombinedList(RhoCandidate* cand, RhoCandList* combinedList, int pdg){
	  /**
	   * @brief: gives back a list of already combined particles
	   * @details: The function creates a list of already combined particles for the analysis
	   */
	  for (int daughter=0; daughter<cand->NDaughters(); daughter++){
		  RhoCandidate * daughterCand = cand->Daughter(daughter);
		  if (daughterCand->PdgCode()==pdg){
			  combinedList->Append(daughterCand);
		  }

	  }

	  combinedList->RemoveClones();


}

void PndTutAnaTask::GetNotCombinedList(RhoCandList combinedList, RhoCandList * candList){
	  for (int j=0; j<combinedList.GetLength(); j++){
		  RhoCandidate * combinedCand = combinedList[j];
		  candList->Remove(combinedCand);
	  }
}


void PndTutAnaTask::numberOfHitsInSubdetector(TString pre, RhoCandidate *c, RhoTuple *n){

	/* This method saves the number of Hits in the MVD, STT and GEM detector
	 * into the RhoTuple.
	 */

	PndPidCandidate *pidCand = (PndPidCandidate*)c->GetRecoCandidate();



	if(pidCand){
		n->Column(pre + "MvdHits", (Int_t) pidCand->GetMvdHits(), 0);
		n->Column(pre + "SttHits", (Int_t) pidCand->GetSttHits(), 0);
		n->Column(pre + "GemHits", (Int_t) pidCand->GetGemHits(), 0);

	}
	else{
		n->Column(pre + "MvdHits", (Int_t) -999, 0);
		n->Column(pre + "SttHits", (Int_t) -999, 0);
		n->Column(pre + "GemHits", (Int_t) -999, 0);
	}
}


void PndTutAnaTask::tagNHits(TString pre, RhoCandidate *c, RhoTuple *n){

	/**@brief Tag the particle with different integers
	 * @details Tag the particle with different integers:
	 * 0: if there is no hit in the detector
	 * 1: sttHits>3 or mvdHits>3 or gemHit>3
	 */

	int tag=0;

	PndPidCandidate * pidCand = (PndPidCandidate*)c->GetRecoCandidate();

	int branch = trackBranch(c);

	if(pidCand){
		int mvdHits = pidCand->GetMvdHits();
		int sttHits = pidCand->GetSttHits();
		int gemHits = pidCand->GetGemHits();


		if(mvdHits>3 || gemHits>3) tag=1;
		else if (sttHits>3 && branch==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack")) tag=1;
		else tag=0;
	}

	n->Column(pre + "HitTag", (Int_t) tag, 0);
}

int PndTutAnaTask::tagHits(RhoCandidate *c){

	/**@brief Tag the particle with different integers
	 * @details Tag the particle with different integers:
	 * 0: if there is no hit in the detector
	 * 1: sttHits>3 or mvdHits>3 or gemHit>3
	 */
	int tag = 0;

	PndPidCandidate * pidCand = (PndPidCandidate*)c->GetRecoCandidate();

	int branch = trackBranch(c);

	if(pidCand){
		int mvdHits = pidCand->GetMvdHits();
		int sttHits = pidCand->GetSttHits();
		int gemHits = pidCand->GetGemHits();

		if(mvdHits>3 || gemHits>3) tag=1;
		else if (sttHits>3 && branch==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack")) tag=1;
		else tag=0;

	}

	return tag;
}

int PndTutAnaTask::trackBranch(RhoCandidate *c){

	int branch=0;

	PndPidCandidate * pid = (PndPidCandidate*)c->GetRecoCandidate();
	if(pid){
		branch = pid->GetTrackBranch();
	}

	return branch;
}

void PndTutAnaTask::TagTrackBranch(RhoCandidate *d0, RhoCandidate *d1, RhoTuple *n){
	/* @brief check if daughter particles cause no hit in the FTS
	 * @details check if daughter particles cause no hit in the FTS. 0 means cause a hit in FTS, 1 means cause no hit in FTS
	 */

	int tagbranch=0;

	PndPidCandidate * pidd0 = (PndPidCandidate*) d0->GetRecoCandidate();
	PndPidCandidate * pidd1 = (PndPidCandidate*) d1->GetRecoCandidate();

	if(pidd0 && pidd1){
		int branchd0 = pidd0->GetTrackBranch();
		int branchd1 = pidd1->GetTrackBranch();

		if(branchd0==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack") & branchd1==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack")){
				tagbranch=1;
		}
	}
	 n->Column("NoFTSHit", (Int_t) tagbranch, -999);
}

void PndTutAnaTask::TagTrackBranch(RhoCandidate *d0, RhoCandidate *d1, RhoCandidate *d2, RhoTuple *n){
	/* @brief check if daughter particles cause no hit in the FTS
	 * @details check if daughter particles cause no hit in the FTS. 0 means cause a hit in FTS, 1 means cause no hin in FTS
	 */

	int tagbranch=0;

	PndPidCandidate * pidd0 = (PndPidCandidate*) d0->GetRecoCandidate();
	PndPidCandidate * pidd1 = (PndPidCandidate*) d1->GetRecoCandidate();
	PndPidCandidate * pidd2 = (PndPidCandidate*) d2->GetRecoCandidate();

	if(pidd0 && pidd1 && pidd2){
		int branchd0 = pidd0->GetTrackBranch();
		int branchd1 = pidd1->GetTrackBranch();
		int branchd2 = pidd2->GetTrackBranch();

		if(branchd0==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack") && branchd1==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack") && branchd2==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack")){
				tagbranch=1;
		}
	}
	 n->Column("NoFTSHit", (Int_t) tagbranch, -999);
}


std::map<int,int> PndTutAnaTask::VertexQaIndex(RhoCandList* candList, float probLimit=0.01){
	  /** @brief  give back the order of the best chi2
	   * @details give back the order of the best chi2!  1 means best, 2: second best (same with negative values for bad chi2 )
	   */

	  std::map<double, int> chi2_good, chi2_bad;

	  for (int j=0; j<candList->GetLength(); ++j){

		  PndKinVtxFitter vtxfitter(candList->Get(j));
		  vtxfitter.Fit();

		  bool failedchi2 = TMath::IsNaN(vtxfitter.GetChi2());
		  bool failedprob = TMath::IsNaN(vtxfitter.GetProb());

		  if(!failedchi2 && !failedprob){

			  if (vtxfitter.GetProb() > probLimit){ //Prob > 0.01
				  chi2_good[vtxfitter.GetChi2()]=j;
			  }
			  else{ //Prob <= 0.01
				  chi2_bad[vtxfitter.GetChi2()]=j;
			  }

		  }
	  }

	  std::map<double, int>::iterator is_good, is_bad;
	  std::map<int, int> indexBestFit;
	  int running = 0;

	  for (is_good = chi2_good.begin(); is_good != chi2_good.end(); is_good++, running++){
		   indexBestFit[is_good->second] = running + 1;
	  }

	  running =0;

	  for (is_bad = chi2_bad.begin(); is_bad != chi2_bad.end(); is_bad++, running++){
		  indexBestFit[is_bad->second] = - (running + 1);
	  }


	  return indexBestFit;
}

std::map<int,int> PndTutAnaTask::MassFitQaIndex(RhoCandList* candList, float m0, float probLimit=0.01){
	  /** @brief  give back the order of the best chi2 for MassFit
	   * @details give back the order of the best chi2 for the MassFit!  1 means best, 2: second best (analoge for bad chi2 with negative values)
	   */

	  if(m0==0) std::cout << "Mass is missing for mass fit" << std::endl;

	  std::map<double, int> chi2_good, chi2_bad;

	  for (int i=0; i<candList->GetLength(); i++){

		  PndKinFitter massfitter(candList->Get(i));
		  massfitter.AddMassConstraint(m0);
		  massfitter.Fit();

		  bool failedchi2 = TMath::IsNaN(massfitter.GetChi2());
		  bool failedprob = TMath::IsNaN(massfitter.GetProb());

		  if(!failedchi2 && !failedprob){

			  if (massfitter.GetProb() > probLimit){
				  chi2_good[massfitter.GetChi2()]=i;
			  }
			  else{
				  chi2_bad[massfitter.GetChi2()]=i;
			  }
		  }
	  }

	  std::map<double,int>::iterator is_good, is_bad;
	  std::map<int,int> bestMassFit;

	  int run =0;

	  for (is_good = chi2_good.begin(); is_good != chi2_good.end(); is_good++, run++){
		  bestMassFit[is_good->second] = run + 1;
	  }

	  run = 0;

	  for (is_bad = chi2_bad.begin(); is_bad != chi2_bad.end(); is_bad++, run++){
		  bestMassFit[is_bad->second] = - (run + 1);
	  }


	  return bestMassFit;
}

void PndTutAnaTask::qaVtxDiff(TString pre, RhoCandidate * c, RhoTuple * n){



	  RhoCandidate * mct = c->GetMcTruth();

	  if(mct){
		  TVector3 v = c->DecayVtx();
		  TVector3 mcv = mct->Daughter(0)->Pos();
		  TVector3 vdiff = v-mcv;
		  TMatrixD cov7 = c->Cov7();

		  n->Column(pre + "rec_truth", 1.0 );
		  n->Column(pre + "diffvx", (Float_t) vdiff.X(), (Float_t)999.);
		  n->Column(pre + "diffvy", (Float_t) vdiff.Y(), (Float_t)999.);
		  n->Column(pre + "diffvz", (Float_t) vdiff.Z(), (Float_t)999.);

		  n->Column(pre + "pullvx", (Float_t) (vdiff.X()/TMath::Sqrt(cov7(0,0))),(Float_t)999.);
		  n->Column(pre + "pullvy", (Float_t) (vdiff.Y()/TMath::Sqrt(cov7(1,1))),(Float_t)999.);
		  n->Column(pre + "pullvz", (Float_t) (vdiff.Z()/TMath::Sqrt(cov7(2,2))),(Float_t)999.);

	  }
	  else{
                  n->Column(pre + "rec_truth", 0.0 );
//		  n->Column(pre + "diffvx", (Float_t) -999.0, 0.0f );
//		  n->Column(pre + "diffvy", (Float_t) -999.0, 0.0f );
//		  n->Column(pre + "diffvz", (Float_t) -999.0, 0.0f );
//
//		  n->Column(pre + "pullvx", (Float_t) -999.0, 0.0f);
//		  n->Column(pre + "pullvy", (Float_t) -999.0, 0.0f);
//		  n->Column(pre + "pullvz", (Float_t) -999.0, 0.0f);
//
	  }
              
}

void PndTutAnaTask::qaMomRes(TString pre, RhoCandidate * c, RhoTuple * n){

	  if(n==0 || c==0) return;

	  RhoCandidate * mct = c->GetMcTruth();
	  float momres = -999.0;

	  if(mct){
		  float p = c->P();
		  float mcp = mct->P();

		  momres = (p-mcp)/mcp;
	  }

	  n->Column(pre + "mom_res", (Float_t) momres, 0.0f);

}

void PndTutAnaTask::col_fit(TString pre, double chi2, double prob, RhoTuple * n){

 n->Column(pre + "chi2", (Float_t) chi2, (Float_t)999.);
 n->Column(pre + "prob", (Float_t) prob, (Float_t)999.);

}

void PndTutAnaTask::col_vtx(TString pre, RhoCandidate * c, RhoTuple * n){
// if(flag==0){
//
// n->Column(pre + "x", -999.0);
// n->Column(pre + "y", -999.0);
// n->Column(pre + "z", -999.0);
//
// }else{ 
 TVector3 Vtx=c->Pos();
 n->Column(pre + "x", (Float_t) Vtx.X(),(Float_t)999. );
 n->Column(pre + "y", (Float_t) Vtx.Y(),(Float_t)999. );
 n->Column(pre + "z", (Float_t) Vtx.Z(),(Float_t)999. );
//}

}





InitStatus PndTutAnaTask::Init(){



	fEvtCount = 0;
  //*** create tuples
  fntpMc = new RhoTuple("ntpMC", "MCTruth info");
  main = new RhoTuple("main", "main info");

  //Create output file for histograms
  //TString outpath;
  //outpath.Append(fOutPath);
  //foutput = TFile::Open( outpath + "output_ana_XiXieta.root","RECREATE");


  // data reader Object
  fAnalysis = new PndAnalysis();
  

  //***Mass selector
  fm0_lambda= TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
  cout<<"Mass of Lambda0: "<<fm0_lambda<<endl;
  lambdaMassSel = new RhoMassParticleSelector("lambda0", fm0_lambda, 0.3);

   fm0_eta = TDatabasePDG::Instance()->GetParticle("eta")->Mass();   // Get nominal PDG mass of the pi0 
   etaMassSel=new RhoMassParticleSelector("eta",fm0_eta,.2);

  fm0_Xi = TDatabasePDG::Instance()->GetParticle("Xi-")->Mass();
  cout<<"Mass of Xi-: "<<fm0_Xi<<endl;
  xiMassSel = new RhoMassParticleSelector("Xi-", fm0_Xi, 0.3);

  fm0_Xis = 2.030;
  cout<<"Mass of Xi(2030)-: "<<fm0_Xis<<endl;
  xi2030MassSel = new RhoMassParticleSelector("Xi(2030)-", fm0_Xis, 0.3);


  TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",3.41706,kFALSE,0.1,0,"",88888);

  fm0_beam = TDatabasePDG::Instance()->GetParticle("pbarpSystem")->Mass();
  cout<<"Mass pbar p system: "<<fm0_beam<<endl;
  

  //*** lorentz vector of the initial particle
//  cout << "BeamMomentum: " << fmom << endl;
  fmom=5.4;
  double p_m0 = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
  fini.SetXYZT(0,0, fmom, sqrt(p_m0*p_m0+ fmom*fmom)+p_m0);
  
//  if(fnevts==0) fnevts = fAnalysis->GetEntries();


  return kSUCCESS;

}

void PndTutAnaTask::Exec(Option_t* op)
{



    int	kPip = 211, kPim = -211,
	kPp = 2212, kaPm = -2212,
	kl0 = 3122, kal0 = -3122,
	kKm = -321,
	kXim = 23314, kaXip = -3312;




	//Dummy RhoCandidate
	RhoCandidate * dummyCand = new RhoCandidate();

	//****************************Analysis*************************************************
//	TVector3 beamBoost = fini.BoostVector();
	PndRhoTupleQA qa(fAnalysis, fmom);
        int i=0,j=0;


         std::map<double, double> chi2map;

        fAnalysis->GetEventInTask();


	 //RhoCandLists for analysis
        RhoCandList protonbar,proton, piplus, piminus,rawgambase,rawgam;
        RhoCandList eta,eta2,lambda0,lambda0bar,ksim,ksip,ppsys;
        RhoCandList vf_lambda0,vf_lambda0bar,vf_eta, eta_vtxfit_list,vf_ksip,vf_ksim,vf_ksim2030;
        RhoCandList reco_pi0,reco_lambda0,reco_lambda0bar,reco_eta,reco_eta2,reco_ksim,reco_ksip,reco_ksim2030,reco_ksip2030,ppsys_4cfit,ppsys_tf,XiXisys,ksi2030_raw,ppsys_raw,ksim_tf,ksip_tf,eta_tf,ppsys_vf,ksim_4cfit,ksip_4cfit,eta_4cfit;
        RhoCandList pip_eta,pim_eta,raw_eta,raw_ksim,raw_ksip,raw_ksim2030,piplus_sur,piminus_sur,piplus_sur2,piminus_sur2,raw_ppsys;


	RhoCandList mclist, all;


        if (!(++fEvtCount)) cout << "evt "<<fEvtCount<<endl;




//	while (fAnalysis->GetEvent() && ++fEvtCount<fnevts){
		//if ((fEvtCount%100)==0) cout << "evt "<< fEvtCount <<endl;
		 cout << "evt "<< fEvtCount <<endl;
           
                //***get MC list and store info, filter here********//
                int flag_Xi_level1=0;
                int flag_Xib_level1=0;
                int flag_lambda_level2=0;
                int flag_lambdab_level2=0;
                int flag_signalMC=0;

                fAnalysis->FillList(mclist, "McTruth");
               for (j=0;j<mclist.GetLength();++j)
              {
                    RhoCandidate *truth=mclist[j];
                    if(truth->PdgCode()==3312){
                      TLorentzVector P4 = truth->P4();
                      TVector3 p3=P4.Vect();
                      
                      
                      int Nd= truth->NDaughters();

                      if (Nd==2&&((truth->Daughter(0)->PdgCode()==3122&&truth->Daughter(1)->PdgCode()==-211)||(truth->Daughter(1)->PdgCode()==3122&&truth->Daughter(0)->PdgCode()==-211))) {
                       flag_Xi_level1=1;
                       TVector3 pos=truth->Daughter(1)->Pos();
                       fntpMc->Column("Angle_Xi_d1pos",(Float_t)p3.Angle(pos),(Float_t)-999.);
                       }  
                  }



                    if(truth->PdgCode()==-3312){
                      int Nd= truth->NDaughters();
                      if (Nd==2&&((truth->Daughter(0)->PdgCode()==-3122&&truth->Daughter(1)->PdgCode()==211)||(truth->Daughter(1)->PdgCode()==-3122&&truth->Daughter(0)->PdgCode()==211))) flag_Xib_level1=1;
                  }


                    if(truth->PdgCode()==3122){
                      int Nd= truth->NDaughters();
                      if (Nd==2&&((truth->Daughter(0)->PdgCode()==2212&&truth->Daughter(1)->PdgCode()==-211)||(truth->Daughter(1)->PdgCode()==2212&&truth->Daughter(0)->PdgCode()==-211))) flag_lambda_level2=1;
                  }
                    if(truth->PdgCode()==-3122){
                      int Nd= truth->NDaughters();
                      if (Nd==2&&((truth->Daughter(0)->PdgCode()==-2212&&truth->Daughter(1)->PdgCode()==211)||(truth->Daughter(1)->PdgCode()==-2212&&truth->Daughter(0)->PdgCode()==211))) flag_lambdab_level2=1;
                  }
                }

               if(flag_Xi_level1==1&&flag_Xib_level1==1&&flag_lambda_level2==1&&flag_lambdab_level2==1) flag_signalMC=1;
                fntpMc->Column("flag_signal", (Int_t)flag_signalMC);
                fntpMc->DumpData();



		//***Setup event shape object

               TString PidSelection = "PidAlgoIdealCharged";//"PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc;PidAlgoEmcBayes";//

		fAnalysis->FillList(all, "All", PidSelection);
		PndEventShape evsh(all, fini, 0.05, 0.1);

		//***Selection


                fAnalysis->FillList(proton,  "ProtonBestPlus","PidAlgoIdealCharged");//????
                fAnalysis->FillList(protonbar, "ProtonBestMinus","PidAlgoIdealCharged");//????        
                fAnalysis->FillList(piplus,  "PionBestPlus","PidAlgoIdealCharged");
                fAnalysis->FillList(piminus, "PionBestMinus","PidAlgoIdealCharged");
                fAnalysis->FillList(rawgambase,"Neutral","PidAlgoIdealNeutral");

                //if(!(proton.GetLength()>=1&&protonbar.GetLength()>=1&&piplus.GetLength()>=2&&piminus.GetLength()>=2&&rawgambase.GetLength()>=2)) continue;
                if(flag_signalMC==1&&(proton.GetLength()>=1&&protonbar.GetLength()>=1&&piplus.GetLength()>=2&&piminus.GetLength()>=2&&rawgambase.GetLength()>=2)) { 


//	    for (int pip=0; pip<piplus.GetLength(); ++pip){
//	        fntpPiPlus->Column("ev",     (Float_t) fEvtCount);
//	        fntpPiPlus->Column("cand",    (Float_t) pip);
//	        fntpPiPlus->Column("ncand",   (Float_t) piplus.GetLength());
//	        fntpPiPlus->Column("McTruthMatch", (bool) fAnalysis->McTruthMatch(piplus[pip]));
//
//	        fntpPiPlus->Column("TrackBranch", (Int_t) trackBranch(piplus[pip]));
//
//	        qa.qaP4("piplus_", piplus[pip]->P4(), fntpPiPlus);
//	        qa.qaCand("piplus_", piplus[pip], fntpPiPlus);
//
//	        numberOfHitsInSubdetector("piplus_", piplus[pip], fntpPiPlus);
//	        tagNHits("piplus_", piplus[pip], fntpPiPlus);
//
//			RhoCandidate * mother_pip;
//			RhoCandidate * truth = piplus[pip]->GetMcTruth();
//			if (truth) mother_pip = truth->TheMother();
//
//			int moth_pip;
//			if (mother_pip==0x0){
//			  moth_pip = 88888;
//			}
//			else
//			  moth_pip = mother_pip->PdgCode();
//
//			fntpPiPlus->Column("Mother", (Float_t) moth_pip);
//
//			qa.qaMcDiff("piplus_", piplus[pip], fntpPiPlus);
//
//			TLorentzVector l;
//			float costheta = -999.;
//			if(truth!=0x0){
//			  l=truth->P4();
//			  costheta = truth->GetMomentum().CosTheta();
//			  qa.qaCand("piplus_MC_", piplus[pip]->GetMcTruth(), fntpPiPlus);
//			}
//			else{
//			  qa.qaCand("piplus_MC_", dummyCand, fntpPiPlus);
//			}
//
//			qa.qaP4("piplus_MC_", l, fntpPiPlus);
//			fntpPiPlus->Column("piplus_MC_CosTheta", (Float_t) costheta);
//
//	        fntpPiPlus->DumpData();
//	    }
//
//	    for (int pim=0; pim<piminus.GetLength(); ++pim){
//	        fntpPiMinus->Column("ev",     (Float_t) fEvtCount);
//	        fntpPiMinus->Column("cand",    (Float_t) pim);
//	        fntpPiMinus->Column("ncand",   (Float_t) piminus.GetLength());
//	        fntpPiMinus->Column("McTruthMatch", (bool) fAnalysis->McTruthMatch(piminus[pim]));
//
//	        int branch = trackBranch(piminus[pim]);
//	        fntpPiMinus->Column("TrackBranch", (Int_t) branch);
//
//
//	        qa.qaP4("piminus_", piminus[pim]->P4(), fntpPiMinus);
//	        qa.qaCand("piminus_", piminus[pim], fntpPiMinus);
//
//	        numberOfHitsInSubdetector("PiMinus_", piminus[pim], fntpPiMinus);
//	        tagNHits("piminus_", piminus[pim], fntpPiMinus);
//
//			RhoCandidate * mother_pim;
//			RhoCandidate * truth = piminus[pim]->GetMcTruth();
//			if (truth) mother_pim = truth->TheMother();
//
//			int moth_pim;
//			if (mother_pim==0x0){
//			  moth_pim = 88888;
//			}
//			else
//			  moth_pim = mother_pim->PdgCode();
//
//			fntpPiMinus->Column("Mother", (Float_t) moth_pim);
//
//			qa.qaMcDiff("piminus_", piminus[pim], fntpPiMinus);
//
//			TLorentzVector l;
//			float costheta = -999.;
//			if(truth!=0x0){
//			  l=truth->P4();
//			  costheta = truth->GetMomentum().CosTheta();
//			  qa.qaCand("piminus_MC_", piminus[pim]->GetMcTruth(), fntpPiMinus);
//			}
//			else{
//			  qa.qaCand("piminus_MC_", dummyCand, fntpPiMinus);
//			}
//
//			qa.qaP4("piminus_MC_", l, fntpPiMinus);
//			fntpPiMinus->Column("piminus_MC_CosTheta", (Float_t) costheta);
//
//	        fntpPiMinus->DumpData();
//
//
//	    }
//
//	    for (int prot=0; prot<proton.GetLength(); ++prot){
//	        fntpProton->Column("ev",     (Float_t) fEvtCount);
//	        fntpProton->Column("cand",    (Float_t) prot);
//	        fntpProton->Column("ncand",   (Float_t) proton.GetLength());
//	        fntpProton->Column("McTruthMatch", (bool) fAnalysis->McTruthMatch(proton[prot]));
//
//	        qa.qaP4("proton_", proton[prot]->P4(), fntpProton);
//	        qa.qaCand("proton_", proton[prot], fntpProton);
//
//	        numberOfHitsInSubdetector("proton_", proton[prot], fntpProton);
//	        tagNHits("proton_", proton[prot], fntpProton);
//
//			RhoCandidate * truth = proton[prot]->GetMcTruth();
//
//			RhoCandidate * mother_prot;
//			if (truth) mother_prot = truth->TheMother();
//			int moth_prot;
//			if (mother_prot==0x0){
//			  moth_prot = 88888;
//			}
//			else
//			  moth_prot = mother_prot->PdgCode();
//
//			fntpProton->Column("Mother", (Float_t) moth_prot);
//
//			TLorentzVector l;
//			float costheta = -999.;
//			if(truth!=0x0){
//			  l=truth->P4();
//			  costheta = truth->GetMomentum().CosTheta();
//			  qa.qaCand("proton_MC_", proton[prot]->GetMcTruth(), fntpProton);
//			}
//			else{
//			  qa.qaCand("proton_MC_", dummyCand, fntpProton);
//			}
//
//			qa.qaP4("proton_MC_", l, fntpProton);
//			fntpProton->Column("proton_MC_CosTheta", (Float_t) costheta);
//
//	        fntpProton->DumpData();
//	    }
//
//	    for (int aProt=0; aProt<protonbar.GetLength(); ++aProt){
//	        fntpprotonbar->Column("ev",     (Float_t) fEvtCount);
//	        fntpprotonbar->Column("cand",    (Float_t) aProt);
//	        fntpprotonbar->Column("ncand",   (Float_t) protonbar.GetLength());
//	        fntpprotonbar->Column("McTruthMatch", (bool) fAnalysis->McTruthMatch(protonbar[aProt]));
//
//	        qa.qaP4("protonbar_", protonbar[aProt]->P4(), fntpprotonbar);
//	        qa.qaCand("protonbar_", protonbar[aProt], fntpprotonbar);
//
//	        numberOfHitsInSubdetector("protonbar_", protonbar[aProt], fntpprotonbar);
//	        tagNHits("protonbar_", protonbar[aProt], fntpprotonbar);
//
//			RhoCandidate * mother_aprot;
//			RhoCandidate * truth = protonbar[aProt]->GetMcTruth();
//			if (truth) mother_aprot = truth->TheMother();
//
//			int moth_aprot;
//			if (mother_aprot==0x0){
//			  moth_aprot = 88888;
//			}
//			else
//			  moth_aprot = mother_aprot->PdgCode();
//
//			fntpprotonbar->Column("Mother", (Float_t) moth_aprot);
//
//
//			TLorentzVector l;
//			float costheta = -999.;
//			if(truth!=0x0){
//			  l=truth->P4();
//			  costheta = truth->GetMomentum().CosTheta();
//			  qa.qaCand("protonbar_MC_", protonbar[aProt]->GetMcTruth(), fntpprotonbar);
//			}
//			else{
//			  qa.qaCand("protonbar_MC_", dummyCand, fntpprotonbar);
//			}
//
//			qa.qaP4("protonbar_MC_", l, fntpprotonbar);
//			fntpprotonbar->Column("protonbar_MC_CosTheta", (Float_t) costheta);
//
//	        fntpprotonbar->DumpData();
//	    }


cout<<__LINE__<<endl;

	    //***Lambda0 -> PiMinus + Proton

	     lambda0.Combine(piminus,proton);
	     lambda0.SetType(kl0);
             chi2map.clear();
                // ***
                // *** do VERTEX FIT (lambda0)
                // ***


                for (j=0;j<lambda0.GetLength();++j)
                {       
                   int flag_l=0;
                   main->Column("lambda0m_raw", (Float_t)lambda0[j]->M(),(Float_t)999.);
                  if (fAnalysis->McTruthMatch(lambda0[0])){
                   flag_l=1;
                   main->Column("mct_lambda0", (Int_t)flag_l);
                    }else{
                   
                  } 



                        PndKinVtxFitter vtxfitter(lambda0[j]);  // instantiate a vertex fitter
                        vtxfitter.Fit();
                        double chi2_vtx = vtxfitter.GetChi2();  // access chi2 of fit
                        double prob_vtx = vtxfitter.GetProb();  // access probability of fit

                        if ( prob_vtx > 0.0 )                          // when good enough, fill some histos
                        {       
                                RhoCandidate *LambdaVfit = lambda0[j]->GetFit();        // access the fitted cand
                                TVector3 Vtx=LambdaVfit->Pos();         // and the decay vertex position
                                col_fit("lambda0_vtf", chi2_vtx,  prob_vtx, main);
                               
 
                                vf_lambda0.Add(LambdaVfit);
                                main->Column("lambda0m_vx", (Float_t)LambdaVfit->M(),(Float_t)999.);
                                int flag2=0;
                                if (fAnalysis->McTruthMatch(LambdaVfit)){
                                 flag2=1;
                                 main->Column("mct_lambda0_vx", (Int_t)flag2);
                                  }               
                        }

                }
                
               vf_lambda0.Select(lambdaMassSel);   

                for (j=0;j<vf_lambda0.GetLength();++j) // !!!!!!!!!!!!   do this fiting with the selcted cand. in vtx fitting
                {
                     PndKinFitter mfitter(vf_lambda0[j]);
                     
                     mfitter.AddMassConstraint(fm0_lambda);
                     mfitter.Fit();
                     double chi2_m = mfitter.GetChi2();
                     if (chi2_m>0)   // exclude unphysical chi2 values  
                     {
                            double prob_mcf = mfitter.GetProb();
                         if (prob_mcf>0.0)
                         {
                          chi2map[chi2_m] = j;
                          col_fit("lambda0_mcf", chi2_m,  prob_mcf, main);
                         }
                     }
                }



cout<<__LINE__<<endl;
                int id_best_proton_forlambda0=0;
                int id_best_piminus_forlambda0=0;
                int flag_lambda0rec=0;
                int lambda0hittag=0;
                TVector3 Vtx_lambda0; 
                if(chi2map.size()>0)
                {
                  flag_lambda0rec=1;
                  reco_lambda0.Put( vf_lambda0[chi2map.begin()->second]->GetFit());
                  id_best_proton_forlambda0 = reco_lambda0[0]->Daughter(0)->Uid();
                  id_best_piminus_forlambda0 = reco_lambda0[0]->Daughter(1)->Uid();
                  main->Column("mcf_lambdam",(Float_t)reco_lambda0[0]->M(),(Float_t)999.);
	          if(tagHits(reco_lambda0[0]->Daughter(0))==1&&tagHits(reco_lambda0[0]->Daughter(1))==1) lambda0hittag=1;
                  col_vtx("lambda0_v", reco_lambda0[0], main);
                  Vtx_lambda0=reco_lambda0[0]->Pos();
                  if (fAnalysis->McTruthMatch(reco_lambda0[0]))
                  {qaVtxDiff("lambda0", reco_lambda0[0], main);}
 
                 }

                chi2map.clear();
	       main->Column("lambda0_HitTag", (Int_t) lambda0hittag);
                   

cout<<__LINE__<<endl;

	    




	     //***AntiLambda0 -> PiPlus + protonbar
	     lambda0bar.Combine(piplus,protonbar);
	     lambda0bar.SetType(kal0);
             chi2map.clear();
                // ***
                // *** do VERTEX FIT (lambda0)
                // ***

                for (j=0;j<lambda0bar.GetLength();++j)
                {       
                        PndKinVtxFitter vtxfitter(lambda0bar[j]);  // instantiate a vertex fitter
                        vtxfitter.Fit();
                        double chi2_vtx = vtxfitter.GetChi2();  // access chi2 of fit
                        double prob_vtx = vtxfitter.GetProb();  // access probability of fit
                        if ( prob_vtx > 0.01 )                          // when good enough, fill some histos
                        {       
                                RhoCandidate *LambdaVfit = lambda0bar[j]->GetFit();        // access the fitted cand
                                TVector3 Vtx=LambdaVfit->Pos();         // and the decay vertex position
                                col_fit("lambda0bar_vtf", chi2_vtx,  prob_vtx, main);
                                vf_lambda0bar.Add(LambdaVfit);
                         
                        }

                }

	        vf_lambda0bar.Select(lambdaMassSel);

                for (j=0;j<vf_lambda0bar.GetLength();++j) // !!!!!!!!!!!!   do this fiting with the selcted cand. in vtx fitting
                {

                     PndKinFitter mfitter(vf_lambda0bar[j]);
                     mfitter.AddMassConstraint(fm0_lambda);
                     mfitter.Fit();
                     double chi2_m = mfitter.GetChi2();
                     if (chi2_m>0)   // exclude unphysical chi2 values  
                     {
                            double prob_mcf = mfitter.GetProb();
                         if (prob_mcf>0.01)
                         {chi2map[chi2_m] = j;
                          col_fit("lambda0bar_mcf", chi2_m,  prob_mcf, main);
                         }
                             

                }
                }    


cout<<__LINE__<<endl;
                int id_best_proton_forlambda0bar=0;
                int id_best_piplus_forlambda0bar=0;
                int flag_lambda0barrec=0;
                int lambda0barhittag=0;
                TVector3 Vtx_lambda0bar;
                if(chi2map.size()>0)
                {
                  flag_lambda0barrec=1;
                  reco_lambda0bar.Put( vf_lambda0bar[chi2map.begin()->second]->GetFit());
                  id_best_proton_forlambda0bar = reco_lambda0bar[0]->Daughter(0)->Uid();
                  id_best_piplus_forlambda0bar = reco_lambda0bar[0]->Daughter(1)->Uid();
	          if(tagHits(reco_lambda0bar[0]->Daughter(0))==1&&tagHits(reco_lambda0bar[0]->Daughter(1))==1) lambda0barhittag=1;
                  col_vtx("lambda0bar_v", reco_lambda0bar[0], main);
                  Vtx_lambda0bar=reco_lambda0bar[0]->Pos();
                  if (fAnalysis->McTruthMatch(reco_lambda0bar[0]))
                  {qaVtxDiff("lambda0bar", reco_lambda0bar[0], main);} 
                  } 
                                  
                chi2map.clear();
                  
	       main->Column("lambda0bar_HitTag", (Int_t) lambda0barhittag);


cout<<__LINE__<<endl;

	      
             //***rec eta->gamma gamma ***********************//
              rawgam.Cleanup();
              double gethr = 0.015;  // GeV 
              for (j=0; j<rawgambase.GetLength(); ++j) {
   
               TLorentzVector rawGamma = rawgambase[j]->P4();
               double e_raw = rawGamma.E();
               double degree_raw = 180*rawGamma.Theta()/TMath::Pi();
               if ( e_raw >= gethr)   rawgam.Add(rawgambase[j]);       // geometry limit on emc forward end-cap
   
               }
   
cout<<__LINE__<<endl;
                 eta.Cleanup();
                 eta.Combine(rawgam,rawgam);
                 eta.SetType(221);
                 eta.Select(etaMassSel);
              int findetapair=0;
              int count_recetaphoton=0;
              for (j=0;j<eta.GetLength();++j)
              {
                  if (fAnalysis->McTruthMatch(eta[j]))
                  { RhoCandidate *mceta_bf = eta[j]->GetMcTruth();
                      if (mceta_bf && mceta_bf->PdgCode()== 221&&mceta_bf->TheMother()->PdgCode()==33314&&findetapair==0)
                      {
                       count_recetaphoton=1;
                       findetapair=1;
                      }
                  }
                 
                  main->Column("count_receta_truth", (Int_t) count_recetaphoton);
  
                  PndKinFitter mfitter_eta(eta[j]);
                  mfitter_eta.AddMassConstraint(fm0_eta);
                  mfitter_eta.Fit();
                  double chi2_eta_m = mfitter_eta.GetChi2();

                 if(chi2_eta_m>0)
                 {
                     double prob = mfitter_eta.GetProb();
                     if (prob>0.0)
                     {
                         chi2map[chi2_eta_m] = j ;
                         col_fit("eta_mcf", chi2_eta_m,  prob, main);
                     }
                  }
                }
            
               int count_receta=0; 
               int count_receta_truth=0; 
              
               TMatrixD p7_eta(7,1);
               TMatrixD cov7_eta(7,7);
               if(chi2map.size()>0)
                {
                  reco_eta.Put(eta[chi2map.begin()->second]->GetFit());

                  p7_eta[0][0]=reco_eta[0]->Px();
                  p7_eta[1][0]=reco_eta[0]->Py();
                  p7_eta[2][0]=reco_eta[0]->Pz();
                  p7_eta[3][0]=reco_eta[0]->E();
                  p7_eta[4][0]=(reco_eta[0]->Pos()).X();
                  p7_eta[5][0]=(reco_eta[0]->Pos()).Y();
                  p7_eta[6][0]=(reco_eta[0]->Pos()).Z();
                  TMatrixD p2Cov=reco_eta[0]->Cov7(); 

                //Change to px,py,pz,E,x,y,z
                 for( i=0; i<7; i++) {
                   for( j=0; j<7; j++) {
                     if(i>=3) {
                       if(j>=3) {
                         cov7_eta[i-3][j-3] = p2Cov[i][j];
                       } else { cov7_eta[i-3][j+4] = p2Cov[i][j]; }
                     } else {
                       if(j>=3) {
                         cov7_eta[i+4][j-3] = p2Cov[i][j];
                       } else { cov7_eta[i+4][j+4] = p2Cov[i][j]; }
                     }
                   }
                 }


                  count_receta=1; 
                  if (fAnalysis->McTruthMatch(reco_eta[0]))
                  { RhoCandidate *mceta = reco_eta[0]->GetMcTruth();
                      if (mceta && mceta->PdgCode()== 221)
                      {
                        
                        count_receta_truth=1;
                       }

                  }
                }
               
               chi2map.clear();
               main->Column("rectag",       (Int_t) count_receta        );
               main->Column("rectag_truth", (Int_t) count_receta_truth  );

               
               
               
               
               


	      //***rec for ksi-***////////////////////////////////////////////////


              piminus_sur.Cleanup();
              ksim.Cleanup();
              vf_ksim.Cleanup();
              raw_ksim.Cleanup();
              reco_ksim.Cleanup();
              chi2map.clear();
  
              //piplus_sur2.Clearup();
  
  
               for (j=0;j<piminus.GetLength();++j)
               {
                 int id_pim = piminus[j]->Uid();
                 if (id_best_piminus_forlambda0!=id_pim)
                 {
                   piminus_sur.Add( piminus[j]);
                  }
  
               }



  
               ksim.Combine(reco_lambda0,piminus_sur);
               ksim.SetType(3312);





               for (j=0;j<ksim.GetLength();++j)
               {
                   PndKinVtxFitter vtxfitter_ksim(ksim[j]);        // instantiate a vertex fitter
                   vtxfitter_ksim.Fit();
                   double chi2_vtx_ksim = vtxfitter_ksim.GetChi2();           // access chi2 of fit
                   if (chi2_vtx_ksim>0)
                   {
                       RhoCandidate *ksiVfit = ksim[j]->GetFit(); 
                       double prob_vtx_ksi = vtxfitter_ksim.GetProb();
                       if (prob_vtx_ksi>0.0 )                  // when good enough, fill some histos
                       {
//                        chi2map[chi2_vtx_ksim] = j;
                        col_fit("ksim_vtf", chi2_vtx_ksim, prob_vtx_ksi, main);
                        vf_ksim.Add(ksiVfit);

                        }
                        }
                        
                    }


                for (j=0;j<vf_ksim.GetLength();++j) // !!!!!!!!!!!!   do this fiting with the selcted cand. in vtx fitting
                {
                     PndKinFitter mfitter(vf_ksim[j]);

                     mfitter.AddMassConstraint(fm0_Xi);
                     mfitter.Fit();
                     double chi2_m = mfitter.GetChi2();
                //     cout<<"mcf Xi chi2"<<chi2_m<<endl;
                     if (chi2_m>0)   // exclude unphysical chi2 values
                     {
                            double prob_mcf = mfitter.GetProb();
                         if (prob_mcf>0.0)
                         {
                          chi2map[chi2_m] = j;
                          col_fit("ksim_mcf", chi2_m,  prob_mcf, main);
                         }
                     }
                }


                int flag_ksimrec=0;
                int ksimhittag=0;                            
                TVector3 Vtx_ksim;
                double ctau_lambda0=999.; 


              if(chi2map.size()>0)
              {

               flag_ksimrec=1;
               reco_ksim.Put(vf_ksim[chi2map.begin()->second]->GetFit());
               main->Column("ksim_m",(Float_t)reco_ksim[0]->M(),(Float_t)999.);
	          if(tagHits(reco_ksim[0]->Daughter(1))==1) ksimhittag=1;
                  col_vtx("ksim_v", reco_ksim[0], main);
                  Vtx_ksim=reco_ksim[0]->Pos();
                  if (fAnalysis->McTruthMatch(reco_ksim[0]))
                  {qaVtxDiff("ksim", reco_ksim[0], main); }
               reco_ksim.Select(xiMassSel);
                }
              
              chi2map.clear();

              if(flag_ksimrec==1&&flag_lambda0rec==1){
              TVector3 Length_lambda0=Vtx_lambda0-Vtx_ksim;
              ctau_lambda0 =Length_lambda0.Mag()*reco_lambda0[0]->M()/reco_lambda0[0]->P(); 
              }

                
              main->Column("lambda0_ctau",(Float_t)ctau_lambda0 );
              main->Column("ksim_HitTag",(Int_t)ksimhittag );


	      //***rec for ksi+***////////////////////////////////////////////////


              ksip.Cleanup();
              piminus_sur.Cleanup();
              vf_ksip.Cleanup();
              raw_ksip.Cleanup();
              reco_ksip.Cleanup();
              chi2map.clear();
  
              for (j=0;j<piplus.GetLength();++j)
                {
                  int id_pip = piplus[j]->Uid();
                  if (id_best_piplus_forlambda0bar!=id_pip)
                  {
                    piplus_sur.Add( piplus[j] );
                   }

                }

               ksip.Combine(reco_lambda0bar,piplus_sur);
               ksip.SetType(-3312);


               for (j=0;j<ksip.GetLength();++j)
               {
                   PndKinVtxFitter vtxfitter_ksip(ksip[j]);        // instantiate a vertex fitter
                   vtxfitter_ksip.Fit();
                   double chi2_vtx_ksip = vtxfitter_ksip.GetChi2();           // access chi2 of fit
                   if (chi2_vtx_ksip>0)
                   {
                       RhoCandidate *ksiVfit = ksip[j]->GetFit();
                       double prob_vtx_ksi = vtxfitter_ksip.GetProb();
                       if (prob_vtx_ksi>0.0 )                  // when good enough, fill some histos
                       {
//                        chi2map[chi2_vtx_ksip] = j;
                        col_fit("ksip_vtf", chi2_vtx_ksip, prob_vtx_ksi, main);
                        vf_ksip.Add(ksiVfit);
                        }
                    }
                        
                }


                for (j=0;j<vf_ksip.GetLength();++j) // !!!!!!!!!!!!   do this fiting with the selcted cand. in vtx fitting
                {
                     PndKinFitter mfitter(vf_ksip[j]);

                     mfitter.AddMassConstraint(fm0_Xi);
                     mfitter.Fit();
                     double chi2_m = mfitter.GetChi2();
                //     cout<<"mcf Xi chi2"<<chi2_m<<endl;
                     if (chi2_m>0)   // exclude unphysical chi2 values
                     {
                            double prob_mcf = mfitter.GetProb();
                         if (prob_mcf>0.0)
                         {
                          chi2map[chi2_m] = j;
                          col_fit("ksip_mcf", chi2_m,  prob_mcf, main);
                         }
                     }
                }


                int flag_ksiprec=0;
                int ksiphittag=0;                            
                TVector3 Vtx_ksip;
                double ctau_lambda0bar=999.; 


              if(chi2map.size()>0)
              {
               flag_ksiprec=1;
               reco_ksip.Put(vf_ksip[chi2map.begin()->second]->GetFit());
               main->Column("ksip_m",(Float_t)reco_ksip[0]->M(),(Float_t)999.);
	          if(tagHits(reco_ksip[0]->Daughter(1))==1) ksiphittag=1;
                  col_vtx("ksip_v", reco_ksip[0], main);
                  Vtx_ksip=reco_ksip[0]->Pos();
                 if (fAnalysis->McTruthMatch(reco_ksip[0]))
                  {qaVtxDiff("ksip", reco_ksip[0], main); }
               reco_ksip.Select(xiMassSel);
                }
              
              chi2map.clear();

              if(flag_ksiprec==1&&flag_lambda0barrec==1){
              TVector3 Length_lambda0bar=Vtx_lambda0bar-Vtx_ksip;
              ctau_lambda0bar=Length_lambda0bar.Mag()*reco_lambda0bar[0]->M()/reco_lambda0bar[0]->P(); 
               }

              main->Column("lambda0bar_ctau",(Float_t)ctau_lambda0bar );
              main->Column("ksip_HitTag",(Int_t)ksiphittag );
               

//	    ******* Vertex of Xi+ Xi- System*****************************
             XiXisys.Cleanup();
             XiXisys.Combine(reco_ksip,reco_ksim);
            TVector3 Vtx_XiXi;
            double ctau_Xim=999.;
            double ctau_Xip=999.;
            TMatrixD Tal0(21,1);
            TMatrixD TV_al0(21,21);
            TV_al0.Zero();     


             for (j=0;j<XiXisys.GetLength();++j)
            {

                  PndKinVtxFitter vtxfitter(XiXisys[j]);        // instantiate a vertex fitter
                   vtxfitter.Fit();
                   double chi2_vtx = vtxfitter.GetChi2();           // access chi2 of fit
                 //   cout<<"XiXi vtf Chi2: "<<chi2_vtx<<endl;
                   if (chi2_vtx>0)
                   {
                       double prob_vtx = vtxfitter.GetProb();
                       if (prob_vtx>0.0 )                  // when good enough, fill some histos
                       {
                        chi2map[chi2_vtx] = j;
                        col_fit("XiXi_vtf", chi2_vtx, prob_vtx, main);
                        Vtx_XiXi=(XiXisys[j]->GetFit())->Pos();

                                TMatrixD Tal0_t=vtxfitter.GetAl0();
              //                  cout<<"Tal0="<<endl;
              //                  Tal0.Print();
                                TMatrixD TV_al0_t=vtxfitter.GetVal0();
              //                  cout<<"TV_al0="<<endl;
              //                  TV_al0.Print();
                             Tal0.SetSub(0,0,Tal0_t);
                             TV_al0.SetSub(0,0, TV_al0_t);
                        }
                    }

                }
             

             
             main->Column("XiXi_vx",(Float_t)Vtx_XiXi.X(),(Float_t)999.);
             main->Column("XiXi_vy",(Float_t)Vtx_XiXi.Y(),(Float_t)999.);
             main->Column("XiXi_vz",(Float_t)Vtx_XiXi.Z(),(Float_t)999.);
 
             chi2map.clear();



//         **********4C fit to*****************************************
                 ppsys_raw.Cleanup();
                 ppsys_raw.Combine(reco_ksip,reco_ksim,reco_eta);// use eta to replace reco_eta
                 ksi2030_raw.Cleanup();
                 ksi2030_raw.Combine(reco_ksim,reco_eta);
                 
                 Tal0.SetSub(14,0,p7_eta);
                 TV_al0.SetSub(14,14,cov7_eta);

                 ppsys_4cfit.Cleanup();
                 ksip_4cfit.Cleanup();
                 ksim_4cfit.Cleanup();
                 eta_4cfit.Cleanup();
                 Int_t mct_ksi2030=0;
                 //TLorentzVector ini(0, 0, 5.2, 6.22224);
                 TLorentzVector ini(0, 0, 5.4, 6.41918);
                for (j=0;j<ppsys_raw.GetLength();++j)
               {
               main->Column("px_t",(Float_t)ppsys_raw[0]->Px(), (Float_t)999.);
               main->Column("py_t",(Float_t)ppsys_raw[0]->Py(), (Float_t)999.);
               main->Column("pz_t",(Float_t)ppsys_raw[0]->Pz(), (Float_t)999.);
               main->Column("E_t",(Float_t)ppsys_raw[0]->E(), (Float_t)999.);
               main->Column("E_t",(Float_t)ppsys_raw[0]->E(), (Float_t)999.);
               main->Column("ksi2030m",(Float_t)ksi2030_raw[0]->M(), (Float_t)999.);
               if(fAnalysis->McTruthMatch(ksi2030_raw[0])) mct_ksi2030=1;
               main->Column("mcflag_ksi2030m",(Int_t)mct_ksi2030);
               cout<<"before 4c"<<endl;
               cout<<"ksim"<<reco_ksim[0]->M()<<endl;
               cout<<"ksip"<<reco_ksip[0]->M()<<endl;
               cout<<"eta"<<reco_eta[0]->M()<<endl;


                 PndKinFitter fitter4c(ppsys_raw[j]);
                 fitter4c.Add4MomConstraint(ini);
                 fitter4c.SetInputMatrix(Tal0, TV_al0);
                 fitter4c.Fit();

                 double chi2 = fitter4c.GetChi2();
                 cout<<"4c chi2="<<chi2<<endl;

                  if (chi2>0)
                  {
                      double prob = fitter4c.GetProb();

                     col_fit("ppsys_4cfit", chi2, prob, main);

                     chi2map[chi2] = j;

                  }

                }
           if(chi2map.size()>0)
           { 

             ppsys_4cfit.Put(ppsys_raw[chi2map.begin()->second]->GetFit());
             ksip_4cfit.Add(ppsys_raw[chi2map.begin()->second]->Daughter(0)->GetFit());
             ksim_4cfit.Add(ppsys_raw[chi2map.begin()->second]->Daughter(1)->GetFit());
             eta_4cfit.Add(ppsys_raw[chi2map.begin()->second]->Daughter(2)->GetFit());


             ///ksip_4cfit.Add(reco_ksip[0]);
             ///ksim_4cfit.Add(reco_ksim[0]);
             ///eta_4cfit.Add(reco_eta[0]);
               cout<<"after 4c"<<endl;
               cout<<"ksim"<<ksim_4cfit[0]->M()<<endl;
               cout<<"ksip"<<ksip_4cfit[0]->M()<<endl;
               cout<<"eta"<<eta_4cfit[0]->M()<<endl;

            main->Column("ppsys_4cfit_mksim", (Float_t)ksim_4cfit[0]->M() , (Float_t)999.);
            main->Column("ppsys_4cfit_mksip", (Float_t)ksip_4cfit[0]->M() , (Float_t)999.);
            main->Column("ppsys_4cfit_eta", (Float_t)eta_4cfit[0]->M() , (Float_t)999.);

            main->Column("ppsys_4cfit_ppsys_px", (Float_t)ppsys_4cfit[0]->Px() , (Float_t)999.);
            main->Column("ppsys_4cfit_ppsys_py", (Float_t)ppsys_4cfit[0]->Py() , (Float_t)999.);
            main->Column("ppsys_4cfit_ppsys_pz", (Float_t)ppsys_4cfit[0]->Pz() , (Float_t)999.);
            main->Column("ppsys_4cfit_ppsys_E", (Float_t)ppsys_4cfit[0]->E() , (Float_t)999.);



              // hetam_mcf_pre_best->Fill(eta_raw[chi2map.begin()->second]->M());

           }

           chi2map.clear();

//~~test/////////////////rec for ksi-(2030)/////////////////////////
           reco_ksim2030.Cleanup();
           reco_ksip2030.Cleanup();
           reco_ksip2030.Combine(ksip_4cfit,eta_4cfit);
           reco_ksim2030.Combine(ksim_4cfit,eta_4cfit);

           reco_ksip2030.SetType(-33314);
           reco_ksim2030.SetType(33314);
           int flag_recoksi2030=0;
           for (j=0;j<reco_ksim2030.GetLength();++j)
          {

           TLorentzVector  p4ksim2030= reco_ksim2030[j]->P4();
           TLorentzVector  p4ksip2030= reco_ksip2030[j]->P4();
           main->Column("ppsys_4cfit_mksimeta", (Float_t)p4ksim2030.M() , (Float_t)999.);
           main->Column("ppsys_4cfit_mksipeta", (Float_t)p4ksip2030.M() , (Float_t)999.);

          if (fAnalysis->McTruthMatch(reco_ksim2030[j])) flag_recoksi2030=1;
           }

          main->Column("count_recoksi2030",(Int_t)flag_recoksi2030  ); 



          main->DumpData();  

/////////////end of fitall()/////////////////////////////////////////



}//if Numbers required

}

void PndTutAnaTask::Finish()
{

	//Write output
	//foutput->cd();
	fntpMc -> GetInternalTree()->Write();
	main->GetInternalTree()->Write();

	//foutput->Save();

	ftimer.Stop();
	Double_t rtime = ftimer.RealTime();
	Double_t ctime = ftimer.CpuTime();
	cout<<endl<<endl;
	cout<<"Macro finisched successfully."<<endl;
	cout<<"Realtime: "<<rtime<<" s, CPU time: "<<ctime<<" s"<<endl;
	cout << "PndTutAnaTask finished" << endl;
}