// *** Mass selector for the J/psi cands
	double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass();  
	RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,0.3);

// *** Mass selector for the pi0 cands
	double m0_pi0 = TDatabasePDG::Instance()->GetParticle("pi0")->Mass(); 
	RhoMassParticleSelector *pi0MassSel = new RhoMassParticleSelector("pi0", m0_pi0, 0.05); 

// *** the lorentz vector of the initial hc
	TLorentzVector ini(0, 0, 5.5827, 6.598952222);

	theAnalysis->FillList(eplus,  "ElectronAllPlus");
	theAnalysis->FillList(eminus, "ElectronAllMinus");
	theAnalysis->FillList(gamma, "Neutral");
		
// *** combinatorics for J/psi -> e+ e-
	jpsi.Combine(eplus, eminus);

// *** combinatorics for pi0 -> gamma gamma
	pi0.Combine(gamma, gamma);

// *** jspi mass selection
        jpsi.Select(jpsiMassSel);
// *** pi0  mass selection
	pi0.Select(pi0MassSel);

// *** combinatorics for hc -> J/psi pi0 pi0
	hc.Combine(jpsi, pi0, pi0);


// ***
// *** do 4C FIT (initial hc system)
// ***
	for (j = 0; j < hc.GetLength(); ++j)
	{
		PndKinFitter fitter(hc[j]);	// instantiate the kin fitter in hc
		fitter.Add4MomConstraint(ini);	// set 4 constraint
		fitter.Fit();		            // do fit

		double chi2_4c = fitter.GetChi2();	// get chi2 of fit
		double prob_4c = fitter.GetProb();	// access probability of fit
		hhc_chi2_4c->Fill(chi2_4c);
		hhc_prob_4c->Fill(prob_4c);
		RhoCandidate *hc_4cf = hc[j]->GetFit(); //get fitted hc itself
		hhc_mass_4cf_prob0->Fill(hc_4cf->M());
        }