#ifndef ___RECO_MESONS___
#define ___RECO_MESONS___

#define TIMEOFFSET 0

//--------------------------------------------------------------
// reconstruction of mother rho,omega,phi meson  (algorithm by Krzysztof Piasecki)
// by using formation times and mother ids of the daughters
// recoverRhoMesons   ( obj ); // has to be  first;
// recoverOmegaMesons ( obj );
// recoverPhiMesons   ( obj );


Int_t  recoverRhoMesons   ( objects& obj ); // have to be  first
Int_t  recoverOmegaMesons ( objects& obj );
Int_t  recoverPhiMesons   ( objects& obj );
//--------------------------------------------------------------


Int_t recoverPhiMesons ( objects& obj )
{
    Int_t N_phis_found = 0;
    smash_p* p1, * p2 ;
    Bool_t found;
    vector<smash_p>& vpart = obj.particles;
    sort(vpart.begin(),vpart.end(),sortFormTime);

    for (UInt_t i = 0; i < vpart.size()-1 ; i++) {
	found = kFALSE;
	p1 = &vpart.at (i);
	p2 = &vpart.at (i+1);

	if ( p1->pdg_mother2 != 0) continue;
	if ( p2->pdg_mother2 != 0) { i++; continue; }

	if ( p1->pdg_mother1 == 333 && p2->pdg_mother1 == 333 && p1->form_time == p2->form_time )
	{
	    if ( (p1->pdg == 321 && p2->pdg == -321) || (p1->pdg == -321 && p2->pdg == 321) ) { found = kTRUE; } // Search for phi -> K+K-
	    if ( (p1->pdg == 311 && p2->pdg == -311) || (p1->pdg == -311 && p2->pdg == 311) ) { found = kTRUE; } // Search for phi -> K0 K0bar
	    if ( (p1->pdg == 221 && p2->pdg ==  22 ) || (p1->pdg ==  22  && p2->pdg == 221) ) { found = kTRUE; } // Search for phi -> eta gamma
	    if ( (p1->pdg == 111 && p2->pdg ==  22 ) || (p1->pdg ==  22  && p2->pdg == 111) ) { found = kTRUE; } // Search for phi -> pi0 gamma
	}

	if (found == kTRUE) {
	    N_phis_found ++ ;

	    smash_p phi;
            phi.reset();

	    phi.v = p1->v + p2->v;                   // Note: phi decay time is stored in FormTime field

	    phi.E      = phi.v.E();
	    phi.mass   = phi.v.M();
	    phi.p0     = phi.v.E();
	    phi.px     = phi.v.Px();
	    phi.py     = phi.v.Py();
	    phi.pz     = phi.v.Pz();
	    phi.pdg    = p1->pdg_mother1;

	    phi.geantID   = 55;
	    phi.charge    = 0;
            phi.form_time = p1->form_time + TIMEOFFSET;
	    phi.weight    = 1.;
            phi.state = 2;
	    phi.v.SetXYZM(phi.px,phi.py,phi.pz,phi.mass);

            p1->state *= -1;
            p2->state *= -1;

	    if(params.verbose_flag >= 1) {
		cout << "\n RECO  ==> EvtNr "<<obj.evt.evtNr<<" Phi from  K+K- , K0 K0bar,eta gamma or pi0 gamma found at pos. [" << i << ',' << i+1 << "]\n" ;
		cout << p1->pdg << '\t' << p2->pdg << '\t' << p1->pdg_mother1 << '\t' << p2->pdg_mother1 << '\t' << p1->form_time << ' ' << p2->form_time << endl;
	    }

	    vpart.erase  ( vpart.begin() + i+1 );
	    vpart.erase  ( vpart.begin() + i   );
	    vpart.insert ( vpart.begin() + i   , phi );
	    i++;

	}
    }

    // Search for phi -> rho pi     (assuming rho has been recovered)

    for (Int_t i = 0; i < vpart.size() ; i++) {
	p1 = &vpart.at (i);

	found = kFALSE;
	if ( isRhoReco (p1->geantID) == kFALSE ) continue;

	Int_t j = i;
        if(j + 1 < vpart.size()) j += 1;
	for (; j >= 0; j--) {                // starting from nearest neighbours, including adjacent 'later one' for safety.
	    if ( j == i ) continue;
	    if (j < 0) break;

	    p2 = &vpart.at (j);

	    if ( p2->pdg_mother1 != 333 )      continue;
	    if ( isPion (p2->pdg) == kFALSE)   continue;
	    if ( p1->charge + p2->charge != 0) continue;        // phi is neutral

	    // We take this pion as the sibling of rho (descending from the same phi)

	    N_phis_found ++ ;

	    smash_p phi;
            phi.reset();

	    phi.v = p1->v + p2->v;                   // Note: phi decay time is stored in FormTime field

	    phi.E      = phi.v.E();
	    phi.mass   = phi.v.M();
            phi.p0     = phi.v.E();
	    phi.px     = phi.v.Px();
	    phi.py     = phi.v.Py();
	    phi.pz     = phi.v.Pz();
	    phi.pdg    = p1->pdg_mother1;

	    phi.geantID   = 55;
	    phi.charge    = 0;
            phi.form_time = p1->form_time + TIMEOFFSET;
	    phi.weight    = 1.;
            phi.state     = 2;
	    phi.v.SetXYZM(phi.px,phi.py,phi.pz,phi.mass);

            p1->state *= -1;
            p2->state *= -1;


	    if(obj.debug >= 1){
		cout << "\n RECO  ==> EvtNr "<<obj.evt.evtNr<<" Phi from  rhopi found at pos. [" << i << ',' << j << "]\n" ;
		cout << p1->pdg << '\t' << p2->pdg << '\t' << p1->pdg_mother1 << '\t' << p2->pdg_mother1 << '\t' << p1->form_time << ' ' << p2->form_time << endl;
	    }

	    if (j == i+1) {
		vpart.erase  ( vpart.begin() + i+1 );
		vpart.erase  ( vpart.begin() + i   );
		vpart.insert ( vpart.begin() + i , phi );
	    } else {
		vpart.erase  ( vpart.begin() + i );
		vpart.erase  ( vpart.begin() + j );
		vpart.insert ( vpart.begin() + j , phi );
	    }

	    break;
	}
    }

    if(N_phis_found>0){
	if(obj.debug >= 2) obj.printEvent(" RECO:  after reconstructing phi");
    }

    return N_phis_found;
}



Int_t recoverOmegaMesons ( objects& obj)
{
    Int_t N_omegas_found = 0;
    smash_p* p1, * p2 ;
    Bool_t found;
    vector<smash_p>& vpart = obj.particles;
    sort(vpart.begin(),vpart.end(),sortFormTime);

    for (UInt_t i = 0; i < vpart.size()-1 ; i++) {
	found = kFALSE;
	p1 = &vpart.at (i);
	p2 = &vpart.at (i+1);

	if ( p1->pdg_mother2 != 0) continue;
	if ( p2->pdg_mother2 != 0) { i++; continue; }

	if ( p1->pdg_mother1 == 223 && p2->pdg_mother1 == 223 && p1->form_time == p2->form_time )
	{
	    if ( (p1->pdg == 111 && p2->pdg == 22) || (p1->pdg == 22 && p2->pdg == 111) ) { found = kTRUE; }// Search for omega -> pi0 gamma
	    if ( isPion (p1->pdg) && isPion (p2->pdg) )                                   { found = kTRUE; }// Search for omega -> pi pi.
	}

	if (found == kTRUE) {
	    N_omegas_found ++ ;
	    smash_p omega;
            omega.reset();

	    omega.v = p1->v + (p2->v);                   // Note: omega decay time is stored in FormTime field

	    omega.E      = omega.v.E();
	    omega.mass   = omega.v.M();
	    omega.p0     = omega.v.E();
	    omega.px     = omega.v.Px();
	    omega.py     = omega.v.Py();
	    omega.pz     = omega.v.Pz();
	    omega.pdg    = p1->pdg_mother1;

	    omega.geantID   = 52;
	    omega.charge    = 0;
            omega.form_time = p1->form_time + TIMEOFFSET;
	    omega.weight    = 1.;
            omega.state     = 2;
	    omega.v.SetXYZM(omega.px,omega.py,omega.pz,omega.mass);

            p1->state *= -1;
            p2->state *= -1;

	    if(obj.debug >= 1) {
		cout << "\n RECO  ==> EvtNr "<<obj.evt.evtNr<<" Omega from  pi0 gamma or pi pi found at pos. [" << i << ',' << i+1 << "]\n" ;
		cout << p1->pdg << '\t' << p2->pdg << '\t' << p1->pdg_mother1 << '\t' << p2->pdg_mother1 << '\t' << p1->form_time << ' ' << p2->form_time << endl;
	    }

	    vpart.erase  ( vpart.begin() + i+1 );
	    vpart.erase  ( vpart.begin() + i   );
	    vpart.insert ( vpart.begin() + i   , omega );
	    i++;
	}
    }

    // Search for omega -> rho pi     (assuming rho has been recovered)
    for (UInt_t i = 0; i < vpart.size() ; i++) {
	p1 = &vpart.at (i);

	found = kFALSE;
	if ( isRhoReco (p1->geantID) == kFALSE ) continue;

	Int_t j = i;
        if(j + 1 < vpart.size()) j += 1;

	for ( ; j >= 0; j--) {		 // starting from nearest neighbours, including adjacent 'later one' for safety.
	    if ( j == i ) continue;
	    if (j < 0) break;

	    p2 = &vpart.at (j);

	    if ( p2->pdg_mother1 != 223 )      continue;        // daughter of omega
	    if ( isPion (p2->pdg) == kFALSE)   continue;        // p2 is pion
	    if ( p1->charge + p2->charge != 0) continue;        // omega is neutral


	    // We take this pion as the sibling of rho (descending from the same omega)

	    N_omegas_found ++ ;

	    smash_p omega;
            omega.reset();

	    omega.v = p2->v + p1->v;                   // Note: omega decay time is stored in FormTime field

	    omega.E      = omega.v.E();
	    omega.mass   = omega.v.M();
	    omega.p0     = omega.v.E();
	    omega.px     = omega.v.Px();
	    omega.py     = omega.v.Py();
	    omega.pz     = omega.v.Pz();
	    omega.pdg    = p1->pdg_mother1;

	    omega.geantID= 52;
	    omega.charge = 0;
	    omega.weight = 1.;
	    omega.state  = 2;
	    omega.v.SetXYZM(omega.px,omega.py,omega.pz,omega.mass);

	    p1->state *= -1;
            p2->state *= -1;

	    if(p2->form_time < p1->form_time) omega.form_time = p2->form_time + TIMEOFFSET;
            else                              omega.form_time = p1->form_time + TIMEOFFSET;

	    if(obj.debug >= 1){
		cout << "\n RECO  ==> EvtNr "<<obj.evt.evtNr<<" Omega from rhopi found at pos. [" << i << ',' << j << "]\n" ;
		cout << p1->pdg << '\t' << p2->pdg << '\t' << p1->pdg_mother1 << '\t' << p2->pdg_mother1 << '\t' << p1->form_time << ' ' << p2->form_time << endl;
	    }

	    if (j == i+1) {
		vpart.erase  ( vpart.begin() + i+1 );
		vpart.erase  ( vpart.begin() + i   );
		vpart.insert ( vpart.begin() + i , omega );
	    } else {
		vpart.erase  ( vpart.begin() + i );
		vpart.erase  ( vpart.begin() + j );
		vpart.insert ( vpart.begin() + j , omega );
	    }
	    break;
	}
    }

    if(N_omegas_found>0){
	if(obj.debug >= 2) obj.printEvent(" RECO: after reconstructing omega");
    }

    return N_omegas_found;
}

Int_t recoverRhoMesons ( objects& obj)
{
    vector<smash_p>& vpart = obj.particles;
    sort(vpart.begin(),vpart.end(),sortFormTime);

    Int_t N_rhos_found = 0;
    smash_p* p1, * p2 ;


    for (Int_t i = 0; i < vpart.size()-1 ; i++) {
	p1 = &vpart.at (i);
	p2 = &vpart.at (i+1);

	if ( p1->pdg_mother2 != 0) continue;
	if ( p2->pdg_mother2 != 0) { i++; continue; }

	if ( p1->form_time != p2->form_time ) continue;

	if ( isPion(p1->pdg)         == kFALSE || isPion(p2->pdg)         == kFALSE )  continue;
	if ( isRho (p1->pdg_mother1) == kFALSE || isRho (p2->pdg_mother1) == kFALSE )  continue;

	N_rhos_found ++ ;

	smash_p rho;
	rho.reset();
	rho.v = p1->v + p2->v;                   // Note: rho decay time is stored in FormTime field

	rho.E      = rho.v.E();
	rho.mass   = rho.v.M();
	rho.p0     = rho.v.E();
	rho.px     = rho.v.Px();
	rho.py     = rho.v.Py();
	rho.pz     = rho.v.Pz();
	rho.pdg    = p1->pdg_mother1;
        rho.form_time = p1->form_time + TIMEOFFSET;
	rho.geantID= geantRho (p1->pdg_mother1);
	rho.charge = rhoCharge(p1->pdg_mother1);
        rho.weight = 1.;
        rho.state  = 2;
	rho.v.SetXYZM(rho.px,rho.py,rho.pz,rho.mass);

        p1->state *= -1;
        p2->state *= -1;


	if(obj.debug >= 1) {
	    cout << "\n RECO  ==> EvtNr "<<obj.evt.evtNr<<" Rho found at pos. [" << i << ',' << i+1 << "]\n" ;
	    cout << p1->pdg << '\t' << p2->pdg << '\t' << p1->pdg_mother1 << '\t' << p2->pdg_mother1 << '\t' << p1->form_time << ' ' << p2->form_time << endl;
	}

	vpart.erase  ( vpart.begin() + i+1 );
	vpart.erase  ( vpart.begin() + i   );
	vpart.insert ( vpart.begin() + i   , rho );
	i++;
    }

    if(N_rhos_found>0){
	if(obj.debug >= 2) obj.printEvent(" RECO:  after reconstructing rho");
    }

    return N_rhos_found;
}


#endif