#include "hparticletool.h"
#include "hcategorymanager.h"


#include "hlinearcategory.h"


#include "hgeantkine.h"
#include "hgeanttof.h"
#include "hgeantshower.h"
#include "hgeantrpc.h"



#include  "TMath.h"
#include  "TAxis.h"
#include  "TVector3.h"

//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////
// HParticleTool
//
// library of static functions
// look here for auxiliary functions for:
// physics analysis, simulation analysis, ...
///////////////////////////////////////////////////////////

ClassImp(HParticleTool)


HParticleTool::HParticleTool()
{
}

HParticleTool::~HParticleTool()
{

}
Float_t HParticleTool::phiSecToLabDeg(Int_t sec, Float_t phiRad)
{
    // Convert phi (rad) angle from sec to Lab phi in deg

    phiRad *= TMath::RadToDeg();

    switch(sec)
    {
    case 0:
	break;

    case 1:
    case 2:
    case 3:
    case 4:
	phiRad += 60.0f * sec;
	break;

    default:
	phiRad -= 60.0f;
	break;
    }

    return phiRad;
}


Float_t HParticleTool::thetaToLabDeg(Float_t thetaRad)
{
    // Convert theta angle (rad) to the coordinate lab system

    return  TMath::RadToDeg() * thetaRad;
}

Float_t HParticleTool::getOpeningAngle(Float_t phi1,Float_t theta1,Float_t phi2,Float_t theta2)
{
    // phi and theta angles of particle 1 and 2 in lab coordinates [deg]
    // returns opening angle [deg]

  //Temporary vector objects for scalar product computation
  TVector3 vec1;
  TVector3 vec2;

  // convert angles to radians
  phi1   *= TMath::DegToRad() ;
  theta1 *= TMath::DegToRad() ;
  phi2   *= TMath::DegToRad() ;
  theta2 *= TMath::DegToRad() ;


  //Above all angles are in degree! We need them in radians!
  vec1.SetMagThetaPhi(1.0,theta1,phi1);
  vec2.SetMagThetaPhi(1.0,theta2,phi2);
  
  return TMath::RadToDeg() * vec1.Angle(vec2);
}

Int_t HParticleTool::findFirstHitInTof(Int_t trackID,Int_t modeTrack)
{
    //-------------------------------------------------
    // find the first track ID entering the TOF
    // Used to suppress the secondaries created in the
    // TOF itself.
    //        0 = realistic (secondaries included)
    //        1 primary particle is stored
    //        2 (default) the first track number entering the tof in SAME SECTOR is stored
    //        3 as 2 but condition on SAME SECTOR && MOD
    //        4 as 2 but SAME SECTOR && MOD && ROD

    Int_t numTrack = trackID;

    if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers

    //------------------------------------------
    // getting categories
    HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE);
    if(!fGeantKineCat){ return trackID; }

    HLinearCategory* fGeantCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE);
    if(!fGeantCat){	return trackID;	}
    //------------------------------------------

    HGeantTof *poldTof, *pnewTof;
    Int_t first = 0;
    HGeantKine* kine = 0;
    kine = HCategoryManager::getObject(kine,fGeantKineCat,numTrack - 1);
    if(kine){
	first=kine->getFirstTofHit();
	if(first != -1){
	    poldTof = (HGeantTof*)fGeantCat->getObject(first);
	    pnewTof = poldTof;
	} else {
	    ::Error("findFirstHitInTof()","No first tof hit!");
	    return numTrack;
	}
    } else {
	::Error("findFirstHitInTof()","Received Zero pointer for kine index = %i !",numTrack);
	return numTrack;
    }
    if(numTrack != poldTof->getTrack()){
	::Error("findFirstHitInTof()","First tof hit not same trackID!");
	return numTrack;
    }

    //--------------------------------------------------------
    // return the track number for
    // the selected option storeFirstTrack

    //--------------------------------------
    // case 0
    if(modeTrack == 0) return numTrack;

    //--------------------------------------
    // case 1
    if(modeTrack == 1)
    {   // return track number of primary particle
	// of the given track
	kine = (HGeantKine*)fGeantKineCat->getObject(numTrack-1);
	kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
	return kine->getTrack();
    }

    //--------------------------------------
    // case 2 and 3
    kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);
    if(kine->getParentTrack() == 0){return numTrack;} // nothing to do

    Int_t s,m,c;
    s = poldTof->getSector();
    // GEANT 15-22 (TOF), 23-26 (TOFINO)
    // numbering reverse !
    m = 21-poldTof->getModule();
    c = 7 -poldTof->getCell();

    first = 0;
    Int_t tempTrack = numTrack;
    while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
    {
	first = kine->getFirstTofHit();
	if(first != -1)
	{  // track is still in TOF

	    // now we have to check if it is in TOF or TOFINO
	    HGeantTof* gtof = (HGeantTof*)fGeantCat->getObject(first);

	    Int_t s1,m1,c1;
	    s1 = m1 = c1 = 0;

	    m1 = 21 - gtof->getModule();
	    if(m1 >= 0)
	    { // inside TOF
		s1 = gtof->getSector();
		c1 = 7-gtof->getCell();

		if(modeTrack == 2 && s == s1)
		{   // case 2 :: check only sector
		    tempTrack = kine->getTrack();
		    pnewTof   = gtof;
		}
		if(modeTrack == 3 && s == s1 && m == m1)
		{   // case 3 :: check only sector,module
		    tempTrack = kine->getTrack();
		    pnewTof   = gtof;
		}
		else if(modeTrack == 4 && s == s1 && m == m1 && c == c1)
		{   // case 4 :: check for same rod
		    tempTrack = kine->getTrack();
		    pnewTof   = gtof;
		}
		else {
		    // track has no TOF hit any longer
		    // which fulfills the condition
		    break;
		}
	    } else {
		// track is in TOFINO
		break;
	    }
	}
	else {
	    // track has no TOF hit any longer,
	    // so the previous object was the one we
	    // searched  for
	    break;
	}
    }
    return tempTrack;
}
Int_t HParticleTool::findFirstHitShowerInTofino(Int_t trackID, Int_t modeTrack)
{
    //-------------------------------------------------
    // Used to suppress the secondaries created in the
    // SHOWER itself.
    //        0 = realistic (secondaries included)
    //        1 primary particle is stored
    //        2 the first track number entering the SHOWER in SAME SECTOR is stored
    //        3 the first track number entering the TOFINO in SAME SECTOR is stored
    //          or the primary track if no TOFINO was found
    //        4 (default) the first track number entering the TOFINO in SAME SECTOR is stored
    //          or the first track in SHOWER if no TOFINO was found


    Int_t numTrack = trackID;

    if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers

    //--------------------------------------
    // case 0  : in = out
    if(modeTrack == 0) { return numTrack; }

    HGeantShower *poldShower;
    Int_t first = 0;
    Int_t parent= 0;
    //------------------------------------------
    // getting categories
    HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE);
    if(!fGeantKineCat){ return trackID; }

    HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE);
    if(!fGeantShowerCat){	return trackID;	}

    HLinearCategory* fGeantTofCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE);
    if(!fGeantTofCat){	return trackID;	}
    //------------------------------------------

    HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);
    if(kine){
	parent = kine->getParentTrack();

	if( parent == 0) {  return numTrack; } // nothing todo

	first = kine->getFirstShowerHit();
	if(first != -1){
	    poldShower = (HGeantShower*)fGeantShowerCat->getObject(first);
	} else {
	    ::Error("findFirstHitShowerInTofino()","No first shower hit!");
	    return numTrack;
	}

    } else {
	::Error("findFirstHitShowerInTofino()","Received Zero pointer for kine!");
	return numTrack;
    }
    if(numTrack != poldShower->getTrack()){
	::Error("findFirstHitShowerInTofino()","First shower hit not same trackID!");
	return numTrack;
    }

    //--------------------------------------------------------
    // return the track number for
    // the selected option modeTrack


    Int_t s = poldShower->getSector();

    //--------------------------------------
    // case 1  : in -> primary
    if(modeTrack == 1)
    {   // return track number of primary particle
	// of the given track

	kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
	return kine->getTrack();
    }

    //--------------------------------------
    // case 2 and 3  : entering SHOWER/TOF

    Int_t tempTrack = numTrack;
    first           = 0;

    //--------------------------------------
    // case 2  : entering SHOWER
    if(modeTrack == 2)
    {
	while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
	{
	    first = kine->getFirstShowerHit();
	    if(first != -1)
	    {
		// track is still in SHOWER

		HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first);
		Int_t s1 = gshower->getSector();

		if(s == s1)
		{   // check only sector
		    tempTrack = kine->getTrack();
		} else {
		    // track has no SHOWER hit any longer
		    // which fulfills the condition
		    break;
		}
	    } else {
		// track has no SHOWER hit any longer,
		// so the previous object was the one we
		// searched for
		break;
	    }
	}
	return tempTrack;
    }
    //--------------------------------------

    //--------------------------------------
    // case 3  : entering TOFINO
    if(modeTrack >= 3)
    {
	Bool_t foundTof  = kFALSE;
	Int_t tempTrack2 = tempTrack;
	do {
	    first = kine->getFirstTofHit();
	    if(first != -1)
	    { // we are in TOF
		HGeantTof* gtof = (HGeantTof*)fGeantTofCat->getObject(first);
		Int_t s1 = gtof->getSector();
		Int_t m  = gtof->getModule();

		if(s == s1 && m > 21 )
		{   // check only sector + TOFINO
		    foundTof  = kTRUE;
		    tempTrack = tempTrack2;
		}
	    }
	    tempTrack2 = kine->getParentTrack();

	} while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0);

	if(foundTof) { tempTrack += 10000; }
	else {
	    kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);


	    if( modeTrack == 3 ){
		// store primaries if no TOFino was found
		kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
		tempTrack = kine->getTrack() + 100000000;
	    } else if (modeTrack == 4){
		//--------------------------------------
		// recover first particle entering SHOWER
		while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
		{
		    first = kine->getFirstShowerHit();
		    if(first != -1)
		    {
			// track is still in SHOWER
			HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first);
			Int_t s1 = gshower->getSector();

			if(s == s1)
			{   // check only sector
			    tempTrack = kine->getTrack();
			} else {
			    // track has no SHOWER hit any longer
			    // which fulfills the condition
			    break;
			}
		    } else {
			// track has no SHOWER hit any longer,
			// so the previous object was the one we
			// searched for
			break;
		    }
		}
		tempTrack += 100000000;
		//--------------------------------------
	    }
	}
    }
    //--------------------------------------

    return tempTrack;
}

Int_t HParticleTool::findFirstHitShowerInRpc(Int_t trackID,Int_t modeTrack)
{
    //-------------------------------------------------
    // Used to suppress the secondaries created in the
    // SHOWER itself.
    //        0 = realistic (secondaries included)
    //        1 primary particle is stored
    //        2 the first track number entering the SHOWER in SAME SECTOR is stored
    //        3 the first track number entering the RPC in SAME SECTOR is stored
    //          or the primary track if no RPC was found
    //        4 (default) the first track number entering the RPC in SAME SECTOR is stored
    //          or the first track in SHOWER if no RPC was found


    Int_t numTrack = trackID;

    if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers

    //--------------------------------------
    // case 0  : in = out
    if(modeTrack == 0) { return numTrack; }

    HGeantShower *poldShower;
    Int_t first = 0;
    Int_t parent= 0;

    //------------------------------------------
    // getting categories
    HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE);
    if(!fGeantKineCat){ return trackID; }

    HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE);
    if(!fGeantShowerCat){	return trackID;	}

    HLinearCategory* fGeantRpcCat = (HLinearCategory*) HCategoryManager::getCategory(catRpcGeantRaw,kFALSE);
    if(!fGeantRpcCat){	return trackID;	}
    //------------------------------------------

    HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);
    if(kine){
	parent = kine->getParentTrack();

	if( parent == 0) {  return numTrack; } // nothing todo

	first = kine->getFirstShowerHit();
	if(first != -1){
	    poldShower = (HGeantShower*)fGeantShowerCat->getObject(first);
	} else {
	    ::Error("findFirstHitShowerInRpc()","No first shower hit!");
	    return numTrack;
	}

    } else {
	::Error("findFirstHitShowerInRpc()","Received Zero pointer for kine!");
	return numTrack;
    }
    if(numTrack != poldShower->getTrack()){
	::Error("findFirstHitShowerInRpc()","First shower hit not same trackID!");
	return numTrack;
    }

    //--------------------------------------------------------
    // return the track number for
    // the selected option modeTrack


    Int_t s = poldShower->getSector();

    //--------------------------------------
    // case 1  : in -> primary
    if(modeTrack == 1)
    {   // return track number of primary particle
	// of the given track

	kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
	return kine->getTrack();
    }

    //--------------------------------------
    // case 2 and 3  : entering SHOWER/RPC

    Int_t tempTrack = numTrack;
    first           = 0;

    //--------------------------------------
    // case 2  : entering SHOWER
    if(modeTrack == 2)
    {
	while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
	{
	    first = kine->getFirstShowerHit();
	    if(first != -1)
	    {
		// track is still in SHOWER

		HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first);
		Int_t s1 = gshower->getSector();

		if(s == s1)
		{   // check only sector
		    tempTrack = kine->getTrack();
		} else {
		    // track has no SHOWER hit any longer
		    // which fulfills the condition
		    break;
		}
	    } else {
		// track has no SHOWER hit any longer,
		// so the previous object was the one we
		// searched for
		break;
	    }
	}
	return tempTrack;
    }
    //--------------------------------------

    //--------------------------------------
    // case 3  : entering RPC
    if(modeTrack >= 3)
    {
	Bool_t foundRPC  = kFALSE;
	Int_t tempTrack2 = tempTrack;
	do {
	    first = kine->getFirstRpcHit();
	    if(first != -1)
	    { // we are in TOF
		HGeantRpc* grpc = (HGeantRpc*)fGeantRpcCat->getObject(first);
		Int_t s1 = grpc->getSector();

		if(s == s1)
		{   // check only sector + RPC
		    foundRPC  = kTRUE;
		    tempTrack = tempTrack2;
		}
	    }
	    tempTrack2 = kine->getParentTrack();

	} while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0);

	if(foundRPC) { tempTrack += 10000; }
	else {
	    kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);


	    if( modeTrack == 3 ){
		// store primaries if no RPC was found
		kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
		tempTrack = kine->getTrack() + 100000000;
	    } else if (modeTrack == 4){
		//--------------------------------------
		// recover first particle entering SHOWER
		while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
		{
		    first = kine->getFirstShowerHit();
		    if(first != -1)
		    {
			// track is still in SHOWER
			HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first);
			Int_t s1 = gshower->getSector();

			if(s == s1)
			{   // check only sector
			    tempTrack = kine->getTrack();
			} else {
			    // track has no SHOWER hit any longer
			    // which fulfills the condition
			    break;
			}
		    } else {
			// track has no SHOWER hit any longer,
			// so the previous object was the one we
			// searched for
			break;
		    }
		}
		tempTrack += 100000000;
		//--------------------------------------
	    }
	}
    }
    //--------------------------------------

    return tempTrack;
}

Float_t HParticleTool::getInterpolatedValue(TH1* h, Float_t xVal, Bool_t warn)
{
    // retrieve content of 1-dim Hist corresponding to x val.
    // The values will be linear interpolated. If warn == kTRUE
    // warning will be emitted if the xval is out of the range
    // of the hist. In case the value is out of range the value
    // of first or last bin will be returned respectively.
    if(!h) return 0;
    TAxis *pA  = h->GetXaxis();
    Int_t binX = pA->FindBin(xVal);

    // We need to start interpolation from the next-lower bin boundary

    if(binX > 1 && binX <= pA->GetNbins())
    {   // only here it make sense to interpolate
	Int_t binXLower = binX - 1;

	Float_t startVal = h->GetBinContent(binXLower);
	Float_t endVal   = h->GetBinContent(binX);

	Float_t startx = pA->GetBinUpEdge(binXLower);
	Float_t slope  = (endVal - startVal)/pA->GetBinWidth(binX);
	Float_t calc   = startVal + slope * (xVal - startx);

	return calc;
    } else { // noting to interpolate....
	if(binX == 1) {
	    return  h->GetBinContent(1);
	} else if(binX < 1) { // lower boundary violation ...return lowest bin
	    if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too low)");
	    return  h->GetBinContent(1);
	} else { // upper boundary violation ...return highest bin
	    if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too high)");
	    return  h->GetBinContent(pA->GetNbins());
	}
    }
}

Stat_t HParticleTool::getValue(TH1* h,Float_t xVal, Float_t yVal, Float_t zVal)
{
    // retrieve value from Histogram (1,2 and 3 dim) corresponding
    // to x val, y val and z val.  If the values are out of range of
    // the histogram axis the lowest/highest bin of the axis will be
    // used. No interpolation performed.

    if(!h) return 0;

    TAxis *pA = h->GetXaxis();
    Int_t  binX, binY, binZ;
    binY = binZ = 0;
    pA   = h->GetXaxis();
    binX = pA->FindBin(xVal);

    if(binX <= 0) { binX = 1; }
    else {
	if(binX > pA->GetNbins()) binX = pA->GetNbins();
    }

    pA   = h ->GetYaxis();
    binY = pA->FindBin(yVal);
    if(binY <= 0) { binY = 1; }
    else {
	if(binY > pA->GetNbins()) binY = pA->GetNbins();
    }

    pA   = h ->GetZaxis();
    binZ = pA->FindBin(zVal);
    if(binZ <= 0) { binZ = 1; }
    else {
	if(binZ > pA->GetNbins()) binZ = pA->GetNbins();
    }

    return h->GetBinContent(h->GetBin(binX, binY, binZ));
}