#include "hparticlepid.h"
#include "hparticletool.h"
#include "hphysicsconstants.h"

//--------category definitions---------
#include "hparticledef.h"
//-------------------------------------

//-------objects-----------------------

#include "hvirtualcand.h"
#include "hparticlecand.h"
#include "hparticlecandsim.h"
//-------------------------------------



//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////
// HParticlePID
//
///////////////////////////////////////////////////////////

ClassImp(HParticlePID)


HParticlePID::HParticlePID(TString beam)
{
    createPIDS();
    nsigma   = 5;
    beamtime = beam;
    rkchi2cut= 1000;
}

HParticlePID::~HParticlePID()
{

}
void HParticlePID::resetPIDS()
{
    mIdToIndex.clear();
    vparticle .clear();
    vcutList  .clear();

}

void HParticlePID::createPIDS()
{
    //---------------------------------------------------------------------------
    resetPIDS();

    // loose cuts
    addPID(2  ,40 ,2000  ,60 ,2000,0);
    addPID(3  ,40 ,2000  ,85 ,2000,0);
    addPID(5  ,40 ,2000  ,60 ,2000,0);
    addPID(6  ,40 ,2000  ,85 ,2000,0);
    addPID(8  ,75 ,2000  ,48 ,2000,0);
    addPID(9  ,84 ,2000  ,70 ,2000,0);
    addPID(11 ,60 ,2000  ,60 ,2000,0);
    addPID(12 ,85 ,2000  ,85 ,2000,0);
    addPID(14 ,291,2000  ,157,4000,0);
    addPID(45 ,60 ,2000  ,60 ,5000,0);
    addPID(46 ,60 ,2000  ,60 ,6000,0);
    addPID(47 ,60 ,2000  ,60 ,4000,0);
    addPID(49 ,60 ,2000  ,60 ,4000,0);
    //---------------------------------------------------------------------------
}

Bool_t HParticlePID::addPID(Int_t id,Float_t min0,Float_t max0,Float_t min1,Float_t max1,UInt_t version)
{
    // add a new momentum cut for PID id (system 0 : min0,max0, system 1 : min1,max1)
    // under the version number version. Since the version corresponds to
    // the index of the cut set list it should consecutive numbers (0,1,2 ....) with
    // already existing versions

    if(id < 0) {
	Warning ("addPID()","id < 0 Skipped!");
	return kFALSE;
    }

    if(version > vcutList.size()){
	Warning ("addPID()","Wrong version = %i, The next new version should be %li",version,vcutList.size());
	return kFALSE;
    }

    if(version < vcutList.size()) { // version already exists

	map<Int_t,momCut>&  cutlist0  =  vcutList[version].sys[0];
	map<Int_t,momCut>&  cutlist1  =  vcutList[version].sys[1];
	if(cutlist0.find(id) == cutlist0.end())
	{  // add new ID
	    momCut cut;
	    cut.setCut(id,min0,max0);
	    cutlist0[id] = cut;
	    cut.setCut(id,min1,max1);
	    cutlist1[id] = cut;
	    if(find(vparticle.begin(),vparticle.end(), id) == vparticle.end()) vparticle.push_back(id);
	    mIdToIndex[id] = vparticle.size()-1;

	} else {
	    cutlist0[id].setCut(id,min0,max0);
	    cutlist1[id].setCut(id,min1,max1);
	}

    } else {  // create new version

	sysList l;
	vcutList.push_back(l);
	map<Int_t,momCut>&  cutlist0  =  vcutList[version].sys[0];
	map<Int_t,momCut>&  cutlist1  =  vcutList[version].sys[1];
	momCut cut;
	cut.setCut(id,min0,max0);
	cutlist0[id] = cut;
	cut.setCut(id,min1,max1);
	cutlist1[id] = cut;
        if(find(vparticle.begin(),vparticle.end(), id) == vparticle.end()) vparticle.push_back(id);
        mIdToIndex[id] = vparticle.size()-1;
    }

    return kTRUE;
}

Int_t HParticlePID::checkPIDS(HParticleCand* pCand,vector<pidvalues>& vValues,UInt_t version)
{
    // checks for all particle hypothesis with the same polarity
    // delta + norm values for beta and MDC dEdx are calculated
    // for each particle id pidvalues are stored to the vector vValues.
    // The index for the id with the smallest norm is returned.
    // If the candidate does not fullfill the input conditions
    // or cannot be assigned to any id -1 is returned.
    // vValues will be cleared before filling.
    // candidates require inner+out MDC fitted segments, RungeKutta fit
    // and a beta measurement.

    vValues.clear();
    if(pCand == 0) return -1;
    if(version >= vcutList.size()) return -1;

    //--------------------------------------------------------------
    // All the possible hypothesis!
    // 5 sigma cut
    if(pCand->getInnerSegmentChi2() > 0 &&
       pCand->getOuterSegmentChi2() > 0 &&
       pCand->getChi2() > 0          &&
       pCand->getChi2() < rkchi2cut  &&
       pCand->getBeta() > 0          &&
       pCand->getSystemUsed() !=-1
      )
    {

	for(UInt_t i=0;i<vparticle.size();i++)
	{
	    Int_t pid =  vparticle[i];
	    Int_t pol =  HPhysicsConstants::charge(pid) > 0 ? 1 : - 1;
	    Int_t sys =  pCand->getSystemUsed();
	    if(pol != pCand->getCharge()) continue;

	    pidvalues val;
	    val.clear();
	    val.id = pid;

	    momCut& cut = vcutList[version].sys[sys][val.id];
	    val.is = HParticleTool::isParticleBeta(val.id,pCand,nsigma, cut.min,cut.max,val.deltaTime,val.sigmaTime,beamtime);
	    if(val.is)
		val.normTime = fabs(val.deltaTime)/val.sigmaTime;
	    else continue;
	    val.is = HParticleTool::isParticledEdx(val.id,pCand,val.deltadEdx,val.sigmadEdx);
	    if(val.is) {
		val.normdEdx = fabs(val.deltadEdx)/val.sigmadEdx;
		// weighted quality value between time and eloss!
		val.norm     = sqrt(val.normTime*val.normTime + val.normdEdx*val.normdEdx*0.6);
	    }
	    else continue;

	    vValues.push_back(val);
	}
	//--------------------------------------------------------------

	//--------------------------------------------------------------
	// Find the best choice for this particle
	Float_t norm  = 9998;
	Int_t ind   = -1;
	for(UInt_t i=0;i<vValues.size();i++) {
	    pidvalues& val = vValues[i];
	    if(val.norm < norm) { norm = val.norm; ind = i; }
	}
	if(ind < 0) return -1;
        else        return ind;
	//--------------------------------------------------------------

    } else return -1;
    //-------------------------------------------------------------------


}

void HParticlePID::print()
{
    cout<<"----------------------------------------------"<<endl;
    cout<<"HParticlePID : " <<endl;
    cout<<"BeamTime  = "<<beamtime.Data()<<endl;
    cout<<"nsigma    = "<<nsigma<<endl;
    cout<<"rkchi2cut = "<<rkchi2cut<<endl;


    for(UInt_t v=0; v < vcutList.size(); v++)
    {

	map<Int_t,momCut>&  cutlist0  =  vcutList[v].sys[0];
	map<Int_t,momCut>&  cutlist1  =  vcutList[v].sys[1];
        cout<<"version = "<<v<<endl;
	for(UInt_t i=0;i<vparticle.size();i++)
	{
            momCut& cut0 = cutlist0[vparticle[i]];
	    momCut& cut1 = cutlist1[vparticle[i]];

	    cout<<setw(2)<<cut0.id
		<<" :  sys0 = ("<<setw(6)<<cut0.min<<","<<setw(6)<<cut0.max<<")"
		<<" :  sys1 = ("<<setw(6)<<cut1.min<<","<<setw(6)<<cut1.max<<")"
		<<endl;
	}
    }
    cout<<"----------------------------------------------"<<endl;

}