#include <TrkNuclei.h>


////////////////////////////////////////////////////////////
/// functions
////////////////////////////////////////////////////////////

/*
 * First charge-vs-beta parameterization provided by wolfgang
 */
Double_t charge_tracker_beta_wolfgang(Double_t *x, Double_t *par){
    
    Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
    Float_t B0[9] = {0 , 5.21108 , 8.53345 , 7.00253 , 8.34315 , 11.9997 , 11.6236 , 2.59707 , 2000 };
    Float_t B1[9] = {0 , -9.6716 , 9.79003 , 38.0199 , 42.6475 , 29.6548 , 43.2637 , 107.077 , 0 };
    Float_t B2[9] = {0 , 5.70349 , -36.0582 , -73.3067 , -69.9404 , -38.9489 , -60.9272 , -152.77 , 0 };
    Float_t B3[9] = {0 , -0.311956 , 21.7421 , 37.4468 , 33.7097 , 17.7431 , 31.2389 , 75.0796 , 0 };
    Float_t x1[9],y1[9];
    Int_t n1 = 9;

    Float_t charge = 1000.;
    Float_t beta    = (Float_t) x[0];
    Float_t dedxtrk = (Float_t) x[1]; dedxtrk = dedxtrk/1.3;

    if( beta>2. || dedxtrk <=0 )return 0.;// it makes no sense to allow beta>2

    Float_t betahelp = pow(beta, 0.8);
    Float_t ym = dedxtrk*betahelp;
    Float_t xb = beta;
    for ( Int_t jj=0; jj<9; jj++ ){
	x1[jj] = B0[jj]+B1[jj]*xb+B2[jj]*xb*xb+B3[jj]*xb*xb*xb;
	y1[jj] = C0[jj]*C0[jj] ;
    }
    TGraph *gr1 = new TGraph(n1,x1,y1);
    TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
    Float_t chelp = spl1->Eval(ym);
    charge = TMath::Sqrt(chelp);

    if(gr1)gr1->Delete();
    if(spl1)spl1->Delete();

    return (Double_t)charge ;

}
/*
 * First charge-vs-rigidity parameterization provided by wolfgang
 * (adpted from charge-vs-deflection parameterization)
 */
Double_t charge_tracker_rigidity_wolfgang(Double_t *x, Double_t *par){

    Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90  };
    Float_t B0[9] = {0 , 0.867439 , 4.04174 , 9.10558 , 15.3956 , 20.6338 , 25.5426 , 32.1287 , 1000  };
    Float_t B1[9] = {0 , -0.128483 , -3.18795 , -3.38137 , -9.61749 , -6.44302 , -8.23092 , -7.30003 , 0 };
    Float_t B2[9] = {0 , 0.764514 , 19.7088 , 38.8391 , 58.1051 , 58.8905 , 55.1826 , 58.9612 , 0 };
    Float_t B3[9] = {0 , 0.00176465 , -9.23937 , -28.0735 , -50.4548 , -54.3154 , -50.3964 , -82.3765 , 0 };
    Float_t B4[9] = {0 , -0.016259 , 1.32851 , 6.68077 , 14.3786 , 16.3259 , 15.835 , 54.1588 , 0 };
    
    
    Float_t yl,yh,m,b,chelp;
    
//     Float_t x1[9],y1[9];
//     Int_t n1 = 9;
    
    Float_t charge  = 1000.;
    Float_t def     = 1./((Float_t)x[0]);
// WM 05-03-2009 calibration valid only for def>0
            def = fabs(def);
    Float_t dedxtrk = (Float_t) x[1]; dedxtrk = dedxtrk/1.3;
    
     Float_t ym = dedxtrk;
     Float_t xb = def;

     for ( Int_t jj=0; jj<8; jj++ ){
	 yl = B0[jj] + B1[jj]*xb + B2[jj]*xb*xb + B3[jj]*xb*xb*xb + B4[jj]*xb*xb*xb*xb;
	 yh = B0[jj+1] + B1[jj+1]*xb + B2[jj+1]*xb*xb + B3[jj+1]*xb*xb*xb + B4[jj+1]*xb*xb*xb*xb;
	 
//	 if ((yl<ym)&&(ym<yh)){
//WM 05-03-2009  dEdx=0 => charge=0
	 if ((yl<=ym)&&(ym<yh)){
	     m = (C0[jj+1]*C0[jj+1] - C0[jj]*C0[jj]) / (yh - yl);
	     b = (C0[jj]*C0[jj]) - m*yl;
	     chelp = m*ym + b;
	     charge= TMath::Sqrt(chelp);
	 }
     }

     return (Double_t)charge ;

}
////////////////////////////////////////////////////////////
/// INITIALIZATIONS & GLOBAL PARAMETERS
////////////////////////////////////////////////////////////

TrkNuclei_parameters * TrkNuclei_parameters::_parameters = 0; 


/**
 * Delete
 */
void TrkNuclei_parameters::Delete(){
    if(charge_vs_beta)charge_vs_beta->Delete();
    if(charge_vs_rig)charge_vs_rig->Delete();
    _parameters=NULL;
}

/**
 * Set external parameters to default values
 */
void TrkNuclei_parameters::SetDefault(){

    cout <<endl << "---------------------------------------------"<<endl;
    cout << "TrkNuclei --> SET DEFAULT PARAMETERS" << endl;
    cout << "---------------------------------------------"<<endl;

    cout << endl << " Charge vs beta "<<endl;
    cout << " >> Wolfgang <<"<<endl;
    charge_vs_beta = new TF2("fch_vs_beta",charge_tracker_beta_wolfgang,0.,1.1,0.,100.);

    cout << endl << " Charge vs rigidity "<<endl;
    charge_vs_rig = new TF2("fch_vs_rig",charge_tracker_rigidity_wolfgang,0.,10000.,0.,100.);
    cout << " >> Wolfgang <<"<<endl;
    cout << "---------------------------------------------"<<endl;

}
/*
 * Set new charge-vs-beta parameterization.
 * Input function should be a 2D function: the first variable is beta, the second one the dedx.
 */
bool TrkNuclei_parameters::SetZ_Beta(TF2* f){
    if(!f)return false;
    cout << "TrkNuclei --> Setting new charge-vs-beta parameterization"<<endl;
    if(charge_vs_beta)charge_vs_beta->Delete();
    charge_vs_beta = f;
    return true;
};
/*
 * Set new charge-vs-beta parameterization.
 * Input function should be a 2D function: the first variable is beta, the second one the dedx.
 */
bool TrkNuclei_parameters::SetZ_Beta( const char* formula ){
    TF2 *f = new TF2("fch_vs_beta",formula,0.,1.1,0.,100.);
    return SetZ_Beta(f);
};
/*
 * Set new charge-vs-rigidity parameterization.
 * Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx.
 */
bool TrkNuclei_parameters::SetZ_Rigidity(TF2* f){
    if(!f)return false;
    cout << "TrkNuclei --> Setting new charge-vs-rigidity parameterization"<<endl;
    if(charge_vs_rig)charge_vs_rig->Delete();
    charge_vs_rig = f;
    return true;
};
/*
 * Set new charge-vs-rigidity parameterization.
 * Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx.
 */
bool TrkNuclei_parameters::SetZ_Rigidity( const char* formula ){
    TF2 *f = new TF2("fch_vs_rig",formula,0.,1.1,0.,100.);
    return SetZ_Rigidity(f);
};

////////////////////////////////////////////////////////////
/// CLASS IMPLEMENTATION
////////////////////////////////////////////////////////////

/**
 * Set event
 */
bool TrkNuclei::Set(TrkTrack* trk){
    if(!trk)return false;
    if(!nutrk)nutrk = new TrkTrack();
    *nutrk = *trk;
//     for(int i=0; i<6; i++){
// 	dedx_x[i]=trk->dedx_x[i];
// 	dedx_y[i]=trk->dedx_y[i];
// 	maxs_x[i]=trk->GetClusterX_MaxStrip(i);
// 	maxs_y[i]=trk->GetClusterY_MaxStrip(i);
//     }    
    return true;
};
bool TrkNuclei::Set(PamTrack* trk){
    if(!trk)return false;
    return Set(trk->GetTrkTrack());
};
bool TrkNuclei::Set(PamLevel2* l2,int ntr){
    if(!l2)return false;
    if(!l2->GetTrkLevel2())return false;
    if(ntr>l2->GetTrkLevel2()->GetNTracks()||ntr<0)return false;
    return Set(l2->GetTrack(ntr)->GetTrkTrack());
};
/**
 * Reset event
 */
void TrkNuclei::Reset(){
    if(nutrk)nutrk->Clear();
//     for(int i=0; i<6; i++){
// 	dedx_x[i]=0.;
// 	dedx_y[i]=0.;
// 	maxs_x[i]=-1;
// 	maxs_y[i]=-1;
//     }
}
/**
 * Method to evaluate the particle charge in each view, as a function of the particle velocity.
 * @param ip plane (0-5)
 * @param iv view (0=x 1=y)
 * @param beta particle velocity
 */
Float_t TrkNuclei::GetZ_Beta(int ip, int iv, float beta){

    if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip,iv)); 

}
/**
 * Method to evaluate the particle charge in each view, as a function of the particle velocity.
 * @param ip plane (0-5)
 * @param beta particle velocity
 */
Float_t TrkNuclei::GetZ_Beta(int ip, float beta){

    if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip)); 

}
/**
 * Method to evaluate the average particle charge, as a function of the particle velocity.
 * @param beta particle velocity
 */
Float_t TrkNuclei::GetZ_Beta(float beta){

    if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX()); 

}
/**
 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
 * @param ip plane (0-5)
 * @param iv view (0=x 1=y)
 * @param rig particle rigidity
 */
Float_t TrkNuclei::GetZ_Rigidity(int ip, int iv, float rig){

    if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip,iv)); 

}
/**
 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
 * @param ip plane (0-5)
 * @param rig particle rigidity
 */
Float_t TrkNuclei::GetZ_Rigidity(int ip, float rig){

    if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip)); 

}
/**
 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
 * @param rig particle rigidity
 */
Float_t TrkNuclei::GetZ_Rigidity(float rig){

    if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
    return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX()); 

}

ClassImp(TrkNuclei);