//_HADES_CLASS_DESCRIPTION 
////////////////////////////////////////////////////////////////////////////
//*-- AUTHOR : J. Markert
////////////////////////////////////////////////////////////////////////////
// HMdcDeDx
//
// This transformation class calculates the "pseudo dEdx" from t2-t1 (time above threshold)
// of all fired drift cells included in one HMdcSeg. The transformation is performed doing a
// normalization of the measured t2-t1 with impact angle and minimum distance to wire
// (MDCCAL2 parametrization) and the algorithm to normalize the single measurements and
// average over all cells included in the segment.
//-------------------------------------------------------------------------
// Calculation of dEdx:
// Float_t calcDeDx(HMdcSeg*,Float_t*,Float_t*,UChar_t*,Float_t*,UChar_t*,
//		    Int_t vers=2,Int_t mod=2, Bool_t useTruncMean=kTRUE, Float_t truncMeanWindow=-99.,Int_t inputselect)
// calculates dedx of a MDC (or 2) segment by doing normalization for
// impact angle and distance from wire for the fired cells of
// the segment. The measurements are sorted in decreasing order,
// the mean and sigma is calculated. Afterwards the truncated mean
// is calculated and returned as dedx. Mean,sigma,number of wires before
// truncating,truncated mean,truncated mean sigma and number of wires
// after truncating can be obtained by the returned pointers.
// Different methods can be selected:
// vers==0 : Input filling from inner HMdcSeg
// vers==1 : Input filling from outer HMdcSeg
// vers==2 : Input filling from inner+outer HMdcSeg (default)
// mod==0  : calc dedx from 1st module in segment
// mod==1  : calc dedx from 2nd module in segment
// mod==2  : calc dedx from whole segment (default)
// useTruncMean: kTRUE (default) apply truncated Mean
// truncMeanWindow (unit: SIGMA RMS of mean TOT): -99 (default) use standard window
// inputselect==0 : fill from segment (default,wires rejected by fit are missing)
// inputselect==1 : fill from cluster
// returns -99 if nothing was calculated
//-------------------------------------------------------------------------
// Settings of the truncated Mean method:
// The truncated mean is calculated according to the set window (in sigma units)
// (setWindow(Float_t window)) arround the mean. For the calculation of the truncated
// mean only drift cells with a dEdx inside the window arround the mean (calculated from
// all drift cells of the segment) will be taken into account.
// With setMinimumWires(Int_t minwires) the algorithm can be forced to keep
// a minimum number of wires. The method will stop removing drift cells from the active
// sample if the minimum number is reached.
////////////////////////////////////////////////////////////////////////////
#include "hmdcdedx2.h"
#include "hpario.h"
#include "hdetpario.h"
#include "hmessagemgr.h"
#include "hparamlist.h"
#include "hmessagemgr.h"
#include "hmdcdetector.h"
#include "hspectrometer.h"
#include "hmatrixcategory.h"
#include "hlinearcategory.h"
#include "hrecevent.h"
#include "hmdcdef.h"
#include "hmdctrackddef.h"
#include "hmdccal1.h"
#include "hmdcseg.h"
#include "hmdchit.h"
#include "hmdcclusinf.h"
#include "hmdcclusfit.h"
#include "hmdcclus.h"
#include "hmdcwirefit.h"
#include "hmdcsizescells.h"
#include "hpidphysicsconstants.h"
#include "hgeantkine.h"

#include "TStyle.h"
#include "TMath.h"
#include "TRandom.h"
#include "TH1.h"
#include <stdlib.h>

ClassImp(HMdcDeDx2)

Float_t HMdcDeDx2::MaxMdcMinDist[4] = {4.,4.,8.,9.};

HMdcDeDx2::HMdcDeDx2(const Char_t* name,const Char_t* title,
		     const Char_t* context)
    : HParCond(name,title,context)
{
    //
    strcpy(detName,"Mdc");
    clear();
    loccal.set(4,0,0,0,0);
    minimumWires  = 5;
    isInitialized = kFALSE;
    setWindow(1.);
    measurements.Set(MAX_WIRES);
    measurements.Reset(-99);
    module = 2;
    method = 2;
    useCalibration = kTRUE;
    debug = kFALSE;
    for(Int_t i = 0; i < 4; i ++){
	MinDistStep[i] = MaxMdcMinDist[i] / (Float_t) N_DIST;
    }
    AngleStep = 5.;

    memset(&parMax[0][0][0][0],0,sizeof(Double_t) * 6 * 4 * N_ANGLE * N_DIST);

}
HMdcDeDx2::~HMdcDeDx2()
{
    // destructor
}
void HMdcDeDx2::clear()
{
    hefr    = 0;
    status  = kFALSE;
    resetInputVersions();
    changed = kFALSE;
}
void HMdcDeDx2::printParam(TString opt)
{
    // prints the parameters of HMdcDeDx2 to the screen.

    cout<<"HMdcDeDx2:"<<endl;
    if(opt.CompareTo("par") == 0 || opt.CompareTo("all") == 0)
    {
	cout<<"par:"<<endl;

	for(Int_t s = 0; s < 6; s ++){
	    for(Int_t m = 0; m < 4; m ++){
		for(Int_t a = 0; a < N_ANGLE; a ++){
		    for(Int_t d = 0; d < N_DIST; d ++){
			for(Int_t i = 0; i < N_PARAM; i ++){
			    printf("s %i m %i angle %2i dist %2i %1.15e %1.15e %1.15e %1.15e\n",
				   s,m,a,d,
				   par[s][m][a][d][0],par[s][m][a][d][1],par[s][m][a][d][2],par[s][m][a][d][3]);
			}
		    }
		}
	    }
	}
    }
    if(opt.CompareTo("parMax") == 0 || opt.CompareTo("all") == 0)
    {
	cout<<"parMax:"<<endl;

	for(Int_t s = 0; s < 6; s ++){
	    for(Int_t m = 0; m < 4; m ++){
		for(Int_t a = 0; a < N_ANGLE; a ++){
		    for(Int_t d = 0; d < N_DIST; d ++){
			printf("s %i m %i angle %2i dist %2i %7.3f\n",
			       s,m,a,d,
			       parMax[s][m][a][d]);
		    }
		}
	    }
	}
    }
    if(opt.CompareTo("shiftpar") == 0 || opt.CompareTo("all") == 0)
    {
	cout<<"shift par:"<<endl;
	for(Int_t s = 0; s < 6; s ++){
	    for(Int_t m = 0; m < 4; m ++){
		for(Int_t a = 0; a < N_ANGLE; a ++){
		    for(Int_t d = 0; d < N_DIST; d ++){
			for(Int_t i = 0; i < N_PARAM;i ++){
			    printf("s %i m %i angle %2i dist %2i %1.15e %1.15e %1.15e %1.15e\n",
				   s,m,a,d,
				   shiftpar[s][m][a][d][0],shiftpar[s][m][a][d][1],shiftpar[s][m][a][d][2],shiftpar[s][m][a][d][3]);
			}
		    }
		}
	    }
	}
    }
    cout<<"window:"<<endl;
    printf("window %7.3f \n",window);
    cout<<"hefr:"<<endl;
    printf("hefr %7.3f \n",hefr);

}
Bool_t HMdcDeDx2::init(HParIo* inp,Int_t* set)
{
    // intitializes the container from an input
    HDetParIo* input = inp->getDetParIo("HCondParIo");
    if (input) return (input->init(this,set));
    return kFALSE;
}
Bool_t HMdcDeDx2::createMaxPar(Bool_t print)
{
    // create the maximum allowed value of ToT to keep
    // the numerical calulations inside the range of double.
    // The procedure loops over all s,m,abin,dbin to set
    // the boundary automatically. This function needs
    // the parameter array for the fit functions to be initialized
    // before. The boundary value are stored inside the container.
    Double_t dEdX;
    Double_t ToT;
    for(Int_t s = 0; s < 6; s ++){
	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t a = 0; a < N_ANGLE; a ++){
		for(Int_t d = 0; d <N_DIST; d ++){
		    Double_t* p = &par[0][m][a][d][0];

		    ToT  = 1.;
		    dEdX = TMath::Power(10.,(TMath::Power((ToT - p[0]) / p[1],(1. / p[2])))) - p[3];
		    while(TMath::Finite(dEdX))
		    {
                        ToT += 1.;
			dEdX = TMath::Power(10.,(TMath::Power((ToT - p[0]) / p[1],(1. / p[2])))) - p[3];
			if(!TMath::Finite(dEdX)) {
			    parMax[s][m][a][d] = ToT - 1;
			    break;
			}
		    }
		    if(print)
		    {
			cout<<"s " <<s
			    <<" m "<<m
			    <<" a "<<a
			    <<" d "<<d
			    <<" ToTmax " <<parMax[s][m][a][d]
			    <<endl;
		    }
		}
	    }
	}
    }

  return kTRUE;
}
Bool_t HMdcDeDx2::initContainer()
{
    if(!isInitialized)
    {

	catcal     = (HCategory*)((HRecEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcCal1);
	if(!catcal) { Error("initContainer()","HMdcCal1 Category not found in current Event!");}

	cathit     = (HCategory*)((HRecEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcHit);
	if(!cathit) { Error("initContainer()","HMdcHit Category not found in current Event!");}

	catclusinf = (HCategory*)((HRecEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcClusInf);
	if(!catclusinf) { Error("initContainer()","HMdcClusInf Category not found in current Event!");}

	catclus     = (HCategory*)((HRecEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcClus);
	if(!catclus) { Error("initContainer()","HMdcClus Category not found in current Event!");}

   	sizescells = HMdcSizesCells::getExObject();
	if(!sizescells)
	{
	    Error("init()","NO HMDCSIZESCELLS CONTAINER IN INPUT!");
	    return kFALSE;
	}
	//if(!sizescells->hasChanged()) sizescells->initContainer();
	isInitialized = kTRUE;
    }
    return kTRUE;
}
Int_t HMdcDeDx2::write(HParIo* output)
{
    // writes the container to an output
    HDetParIo* out = output->getDetParIo("HCondParIo");
    if (out) return out->write(this);
    return -1;
}
void HMdcDeDx2::putParams(HParamList* l)
{
    // Puts all params of HMdcDeDx2 to the parameter list of
    // HParamList (which ist used by the io);
    if (!l) return;
    l->add("par"      ,&par     [0][0][0][0][0],6 * 4 * N_ANGLE * N_DIST * N_PARAM      );
    l->add("parMax"   ,&parMax  [0][0][0][0]   ,6 * 4 * N_ANGLE * N_DIST                );
    l->add("shiftpar" ,&shiftpar[0][0][0][0][0],6 * 4 * N_ANGLE * N_DIST * N_SHIFT_PARAM);
    l->add("window"   ,window);
    l->add("hefr"     ,hefr);


}
Bool_t HMdcDeDx2::getParams(HParamList* l)
{
    if (!l) return kFALSE;
    if (!( l->fill("par"      ,&par     [0][0][0][0][0],6 * 4 * N_ANGLE * N_DIST * N_PARAM      ))) return kFALSE;
    if (!( l->fill("parMax"   ,&parMax  [0][0][0][0]   ,6 * 4 * N_ANGLE * N_DIST                ))) return kFALSE;
    if (!( l->fill("shiftpar" ,&shiftpar[0][0][0][0][0],6 * 4 * N_ANGLE * N_DIST * N_SHIFT_PARAM))) return kFALSE;
    if (!( l->fill("window"   ,&window                                                          ))) return kFALSE;
    if (!( l->fill("hefr"     ,&hefr                                                            ))) return kFALSE;

    return kTRUE;
}

void HMdcDeDx2::sort(Int_t nHits)
{
    // Puts the measurement values into decreasing order.

    Int_t nhit=nHits;
    if(nHits>MAX_WIRES) nhit=MAX_WIRES;

    register Int_t a,b,c;
    Int_t exchange;
    Double_t val;

    for(a = 0; a < nhit - 1; ++ a)
    {
	exchange = 0;
	c        = a;
	val      = measurements[a];

	for(b = a + 1; b < nhit; ++ b)
	{
	    if(measurements[b] > val)   // : < increasing  :> decreasing
	    {
		c = b;
		val = measurements[b];
		exchange = 1;
	    }
	}
	if(exchange)
	{
	    measurements[c] = measurements[a];
	    measurements[a] = val;
	}
    }
}
Float_t HMdcDeDx2::calcDeDx(HMdcSeg* seg[2],
			    Float_t* meanold,Float_t* sigmaold,UChar_t* nwire,
			    Float_t* sigmanew,UChar_t* nwiretrunc,
			    Int_t vers,Int_t mod,
			    Bool_t useTruncMean,
			    Float_t truncMeanWindow,
			    Int_t inputselect)
{
    // calculates dedx of a MDC segment by doing normalization for
    // impact angle and distance from wire for the fired cells of
    // the segment. The measurements are sorted in decreasing order,
    // the mean and sigma is calculated. Afterwards the truncated mean
    // according to the set window (in sigma units) arround the mean
    // is calculated and returned as dedx. Mean,sigma,number of wires before
    // truncating,truncated mean,truncated mean sigma and number of wires
    // after truncating can be obtained by the returned pointers.
    // Different methods can be selected:
    // vers==0 : Input filling from inner HMdcSeg
    // vers==1 : Input filling from outer HMdcSeg
    // vers==2 : Input filling from inner+outer HMdcSeg (default)
    // mod==0  : calc dedx from 1st module in segment
    // mod==1  : calc dedx from 2nd module in segment
    // mod==2  : calc dedx from whole segment (default)
    // useTruncMean: kTRUE (default) apply truncated Mean
    // truncMeanWindow (unit: SIGMA RMS of mean TOT): -99 (default) use standard window
    // inputselect==0 : fill from segment (default, wires rejected by fit are missing)
    // inputselect==1 : fill from cluster
    // returns -99 if nothing was calculated

    //to be shure initialized values are returned
    *meanold   = -1.;
    *sigmaold  = -1.;
    *nwire     =  0;
    *sigmanew  = -1.;
    *nwiretrunc=  0;
    Float_t dedx       = -99.;
    Float_t mean       = -99.;
    Float_t sigma      = -99.;
    UChar_t nWire      =   0;
    UChar_t nWireTrunc =   0;

    if(!catcal ) return mean; // no cal1 data in input !

    if(!isInitialized) {
	Error("HMdcDeDx2::calcDeDx()","Container HMdcDeDx2 has not been initialized!");
	return mean;

    }
    if(vers >= 0 && vers <= 2) method = vers;

    if(seg[0] == 0 && (method == 0 || method == 2))
    {
	Error("HMdcDeDx2::calcDeDx()","ZERO POINTER FOR inner HMdcSeg RECEIVED!");
	return mean;
    }
    if(seg[1] == 0 && (method == 1))
    {
	return mean;
    }

    if(mod >= 0 && mod < 3)
    {
	module = mod;
    }
    else
    {
	Warning("calcDeDx()","Unknown module type! Will use both modules of Segment!");
        module = 2;
    }

    //---------------------------------------------------
    // getting input
    measurements.Reset(-99);

    nWire = fillingInput(seg,inputselect);
    if(debug) {
	cout<<"calcDeDx() from fillingInput(): "
	    <<"nw "      <<((Int_t)nWire)
	    <<" methode "<<method
            <<" module " <<module
	    <<endl;
    }
    if(nWire == 0) return mean;

    *nwire = nWire;

    if(nWire > 0) mean = TMath::Mean(nWire,measurements.GetArray(), NULL);
    *meanold = mean;
    //---------------------------------------------------
    // calculating sigma
    if(nWire >= 1)
    {
        sigma       = TMath::RMS(nWire,measurements.GetArray());
        *sigmaold   = sigma;
        //--------------------truncating measurements outside +- x*sigma
	if (useTruncMean) {
	    nWireTrunc  = select(mean,sigma,nWire,truncMeanWindow);
	} else {
	    sort(nWire);
            nWireTrunc=nWire;
	}
        *nwiretrunc = nWireTrunc;
        //--------------------calculating truncated mean
        dedx        = TMath::Mean(nWireTrunc,measurements.GetArray(), NULL);
        //--------------------calculating new sigma
        sigma       = TMath::RMS(nWireTrunc,measurements.GetArray());
        *sigmanew   = sigma;

    }else{
	Warning("calcDeDx()","nWire <=1 : skipped %i and %i with t2<=-998 !",
		ctskipmod1,ctskipmod1);
    }

    //---------------------------------------------------
    // now calibrate mean value of ToT to dEdX
    Int_t s = 0;
    if(seg[0]) s = seg[0]->getSec();
    dedx = toTTodEdX(s,ref_MOD,ref_ANGLE,ref_DIST,dedx);
    if(dedx > 1000) {
	if(debug){cout<<"calcDeDx() : dEdX>1000 -->"<<dedx<<endl;}
	dedx = 1000.;
    }
    //---------------------------------------------------

    return dedx;
}
UChar_t HMdcDeDx2::fillingInput(HMdcSeg* seg[2],Int_t inputselect)
{
    // Fills array of measurements from HMdcSeg and returns the number of fired wires in Segment
    // The single measurements are normalized by impact angle and distance from wire. Only the first
    // wires < MAX_WIRES are accepted.
    // inputselect==0 : fill from segment (default,wires rejected by fit are missing)
    // inputselect==1 : fill from cluster

    ctskipmod0 = 0;
    ctskipmod1 = 0;

    Float_t t1,t2;
    Double_t alpha,mindist;
    Double_t x1,y1,z1,x2,y2,z2;

    UChar_t nWire = 0;

    if(!catcal ) return nWire; // no cal1 data in input !

    if(inputselect == 1 &&
       (!cathit || !catclusinf || !catclus)
      ) return nWire; // no cluster data in input !

    HMdcClus* clus[2] = {0};
    //----------------------------------------getting input--------------------------------------------------
    Int_t low = 0;
    Int_t up  = 2;

    if     (method == 0) {low = 0;up = 1;}   // inner seg
    else if(method == 1) {low = 1;up = 2;}   // outer seg
    else if(method == 2) {low = 0;up = 2;}   // both  seg

    for(Int_t io = low; io < up; io ++)
    {
        if(seg[io] == 0) continue;

        //---------------------------------------------------------------------------
	if(inputselect == 1)
	{
	    // we have to get the cluster for looping
	    // the drift cells
	    for(Int_t ihit = 0; ihit < 2; ihit ++)
	    {
		Int_t hitind = seg[io]->getHitInd(ihit);
		if(hitind != -1)
		{   // hit and clusinf have the same index!
		    HMdcClusInf* clusinf = (HMdcClusInf*) catclusinf->getObject(hitind);
		    if(clusinf)
		    {
			clus[io] = (HMdcClus*) catclus->getObject(clusinf->getClusIndex());
			if(clus[io] == 0){
			    Error("fillingInput()","Zero pointer for HMdcClus object retrieved!");
			    return nWire;
			}
		    }
		    else Error("fillingInput()","Zero pointer for HMdcClusInd object retrieved!");

		    break;  // found cluster -> exit the loop
		}

	    }
	}
        //---------------------------------------------------------------------------


	for(Int_t l = 0; l < 12; l ++)
	{
	    if(module == 0 && l > 5) continue;  // fill for 1st module only
	    if(module == 1 && l < 6) continue;  // fill for 2st module only

	    Int_t nCell = 0;
            // input selection
	    if(inputselect == 1){nCell = clus[io]->getNCells(l);}
	    else                {nCell = seg [io]->getNCells(l);}

	    loccal[0]   = seg[io]->getSec();
	    Int_t ioseg = seg[io]->getIOSeg();

	    for(Int_t i = 0; i < nCell; i ++)
	    {
		if(ioseg == 0){
		    (l < 6)? loccal[1] = 0 : loccal[1] = 1;
		}else if(ioseg == 1){
		    (l < 6)? loccal[1] = 2 : loccal[1] = 3;
		}
		(l < 6)? loccal[2] = l : loccal[2] = l - 6;

		// input selection
		if(inputselect == 1){loccal[3] = clus[io]->getCell(l,i);}
                else                {loccal[3] = seg [io]->getCell(l,i);}

		calcSegPoints(seg[io],x1,y1,z1,x2,y2,z2);
		(*sizescells)[loccal[0]][loccal[1]][loccal[2]].getImpact(x1,y1,z1,x2,y2,z2,loccal[3],alpha,mindist);

		HMdcCal1* cal1 = (HMdcCal1*)catcal->getObject(loccal);
		if(cal1 != 0)
		{
		    t1  = cal1->getTime1();
		    t2  = cal1->getTime2();

		    if(t2 != -998 && t2 != -999 &&
                       t1 != -998 && t1 != -999 &&
		       TMath::Finite(t1)        &&
                       TMath::Finite(t2)
		      )
		    {   // some times t1 or t2 can be missing (-999,-998)
			// calculating mean value

			if(nWire < MAX_WIRES - 1)
			{
			    if(useCalibration)measurements[nWire] = normalize(loccal[0],loccal[1],alpha,mindist,t2-t1);
			    else              measurements[nWire] = t2-t1;
			    nWire ++;
			}else{
			    if(debug) Warning("fillingInput()","Number of wires in Segment=%i > MAX_WIRES = %i! Skipped!",nWire,MAX_WIRES);
			}
		    } else {
			if(loccal[1] == 0 || loccal[1] == 2) ctskipmod0 ++;
			if(loccal[1] == 1 || loccal[1] == 3) ctskipmod1 ++;}
		}else{
		    Warning("calcDeDx()","ZERO pointer recieved for cal1 object!");
		}
	    }
	}
    }
    return nWire;
}

UChar_t HMdcDeDx2::select(Float_t mean,Float_t sigma,UChar_t nWire,Float_t wind)
{
    // loops over the measurement array with nWire values and puts values to
    // -99 if the measurements are outside the wanted window arround the mean value mean
    // (set with setWindow() (sigma units)).
    // The procedure keeps at least a number of minimum wires set with setMinimumWires().

    UChar_t nWTrunc = nWire;
    UChar_t count   = 0;
    UChar_t minW    = (UChar_t)getMinimumWires();

    if(nWTrunc < minW)
    {
	if(debug){cout<<"select : skipp because nWT<minW "<<endl;}
	return nWTrunc;
    }

    //---------------------------------------------------
    // decide about the window for truncated Mean method
    Float_t tempWindow = getWindow();
    if(wind > 0) tempWindow = wind;
    //---------------------------------------------------

    //--------------------sorting measurements in decreasing order
    sort(nWire);
    if(debug){
	cout<<"---------------------------------------------------"<<endl;
	cout<<"measurements before truncation :"<<endl;
	for(Int_t i = 0; i < nWire; i ++){
	    cout<<i<<" "<<measurements[i]<<endl;
	}
    }

    Bool_t fail_high = kFALSE;
    Bool_t fail_low  = kFALSE;

    while(nWTrunc > minW &&
	  count < nWire )
    {
	//------------------------------------------------------------
	// starting from highest val
        if(!fail_high && fabs(measurements[count] - mean) > (tempWindow * sigma))
        { // if measurement outside the window
	    measurements[count] = -99;
	    nWTrunc --;
        } else  fail_high = kTRUE;
	//------------------------------------------------------------

	//------------------------------------------------------------
	// taking lowest val next
	if(nWTrunc > minW)
	{
	    if(!fail_low && fabs(measurements[nWire - 1 - count] - mean) > (tempWindow * sigma))
	    { // if measurement outside the window
		measurements[nWire - 1 - count] = -99;
		nWTrunc --;
	    }
	    else  fail_low = kTRUE;
	} else  fail_low = kTRUE;
	//------------------------------------------------------------
	count ++;
        if(fail_low && fail_high) break;
    }
    sort(nWire);

    if(debug){
	cout<<"---------------------------------------------------"<<endl;
	cout<<"measurements after truncation :"<<endl;
	for(Int_t i = 0; i < nWTrunc; i ++){
	    cout<<i<<" "<<measurements[i]<<endl;
	}
	cout<<"method "  <<method
	    <<" window " <<tempWindow
	    <<" mean "   <<mean
	    <<" meanN "  <<TMath::Mean(nWTrunc,measurements.GetArray())
	    <<" sig "    <<sigma
	    <<" sigN "   <<TMath::RMS(nWTrunc,measurements.GetArray())
	    <<" w "      <<((Int_t)nWire)
	    <<" trunc w "<<((Int_t)nWTrunc)
	    <<endl;
    }

    return  nWTrunc;

}
void HMdcDeDx2::setFuncPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t* p,Int_t size)
{
    // set the values for the dedx functions by s, m, angle bin, dist bin. Needs pointer
    // to parameter array[N_PARAM]
    if(size == N_PARAM             &&
       abin >= 0 && abin < N_ANGLE &&
       dbin >= 0 && dbin < N_DIST  &&
       s >= 0 && s < 6             &&
       m >= 0 && m < 4){
	for(Int_t i = 0; i < size; i ++) par[s][m][abin][dbin][i] = p[i];
    } else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncPar(Double_t* p)
{
    // set all values for the dedx functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST][N_PARAM]
   for(Int_t s = 0; s < 6;s ++) {
	for(Int_t m = 0; m < 4; m ++) {
	    for(Int_t a = 0; a < N_ANGLE; a ++) {
		for(Int_t d = 0; d < N_DIST; d ++) {
		    for(Int_t i = 0; i < N_PARAM; i ++) {

			Int_t maxS = 4 * N_ANGLE * N_DIST * N_PARAM;
                        Int_t maxM = N_ANGLE * N_DIST * N_PARAM;
                        Int_t maxA = N_DIST * N_PARAM;

			par[s][m][a][d][i] = p[s * maxS + m * maxM + a * maxA + d * N_PARAM + i];
		    }
		}
	    }
	}
    }
}
void HMdcDeDx2::getFuncPar(Double_t* p)
{
    // set all values for the dedx functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST][N_PARAM]
   for(Int_t s = 0;s < 6; s ++) {
	for(Int_t m = 0; m < 4; m ++) {
	    for(Int_t a = 0; a < N_ANGLE; a ++) {
		for(Int_t d = 0; d < N_DIST; d ++) {
		    for(Int_t i = 0; i < N_PARAM; i ++) {

			Int_t maxS = 4 * N_ANGLE * N_DIST * N_PARAM;
                        Int_t maxM = N_ANGLE * N_DIST * N_PARAM;
                        Int_t maxA = N_DIST * N_PARAM;

			p[s * maxS + m * maxM + a * maxA + d * N_PARAM + i] = par[s][m][a][d][i];
		    }
		}
	    }
	}
    }
}
void HMdcDeDx2::setFuncMaxPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t p)
{
    // set the values for the dedx functions by s, m, angle bin, dist bin.
    if(abin >= 0 && abin < N_ANGLE &&
       dbin >= 0 && dbin < N_DIST  &&
       s >= 0 && s < 6             &&
       m >= 0 && m < 4){
	parMax[s][m][abin][dbin] = p;
    } else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncMaxPar(Double_t* p)
{
    // set all values for the dedx functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST]
    for(Int_t s = 0; s < 6; s ++) {
	for(Int_t m = 0; m < 4; m ++) {
	    for(Int_t a = 0; a < N_ANGLE; a ++) {
		for(Int_t d = 0; d < N_DIST;d ++) {
		    Int_t maxS = 4 * N_ANGLE * N_DIST;
		    Int_t maxM = N_ANGLE * N_DIST;
		    Int_t maxA = N_DIST;

		    parMax[s][m][a][d] = p[s * maxS + m * maxM + a * maxA + d];
		}
	    }
	}
    }
}
void HMdcDeDx2::getFuncMaxPar(Double_t* p)
{
    // set all values for the dedx functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST]
   for(Int_t s = 0; s < 6; s ++) {
       for(Int_t m = 0; m < 4; m ++) {
	   for(Int_t a = 0; a < N_ANGLE; a ++) {
	       for(Int_t d = 0; d <N_DIST; d ++) {
		   Int_t maxS = 4 * N_ANGLE * N_DIST;
		   Int_t maxM = N_ANGLE * N_DIST;
		   Int_t maxA = N_DIST;

		   p[s * maxS + m * maxM + a * maxA + d] = parMax[s][m][a][d];
	       }
	   }
       }
   }
}

void HMdcDeDx2::setFuncWidthPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t* p,Int_t size)
{
    // set the values fpr the dedx width functions by s, m, angle bin, dist bin. Needs pointer
    // to parameter array[N_SHIFT_PARAM]

    if(size == N_SHIFT_PARAM      &&
       abin >= 0 && abin < N_ANGLE&&
       dbin >= 0 && dbin < N_DIST &&
       s >= 0 && s < 6            &&
       m >= 0 && m < 4){
	for(Int_t i = 0; i < size; i ++) shiftpar[s][m][abin][dbin][i] = p[i];
    } else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncWidthPar(Double_t* p)
{
    // set all values for the dedx width functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST][N_SHIFT_PARAM]
    for(Int_t s = 0;s < 6; s ++) {
	for(Int_t m = 0; m < 4; m ++) {
	    for(Int_t a = 0; a < N_ANGLE; a ++) {
		for(Int_t d = 0; d < N_DIST; d ++) {
		    for(Int_t i = 0; i < N_SHIFT_PARAM; i ++) {

			Int_t maxS = 4 * N_ANGLE * N_DIST * N_SHIFT_PARAM;
                        Int_t maxM = N_ANGLE * N_DIST * N_SHIFT_PARAM;
                        Int_t maxA = N_DIST * N_SHIFT_PARAM;

			shiftpar[s][m][a][d][i] = p[s * maxS + m * maxM + a * maxA + d * N_SHIFT_PARAM + i];
		    }
		}
	    }
	}
    }
}
void HMdcDeDx2::getFuncWidthPar(Double_t* p)
{
    // set all values for the dedx width functions. Needs pointer
    // to parameter array[6][4][N_ANGLE][N_DIST][N_SHIFT_PARAM]
    for(Int_t s = 0; s < 6; s ++) {
	for(Int_t m = 0; m < 4; m ++) {
	    for(Int_t a = 0; a < N_ANGLE; a ++) {
		for(Int_t d = 0; d < N_DIST; d ++) {
		    for(Int_t i = 0; i < N_SHIFT_PARAM; i ++) {

			Int_t maxS = 4 * N_ANGLE * N_DIST * N_SHIFT_PARAM;
                        Int_t maxM = N_ANGLE * N_DIST * N_SHIFT_PARAM;
                        Int_t maxA = N_DIST * N_SHIFT_PARAM;

			p[s * maxS + m * maxM + a * maxA + d * N_SHIFT_PARAM + i] = shiftpar[s][m][a][d][i];
		    }
		}
	    }
	}
    }
}


void  HMdcDeDx2::calcSegPoints(HMdcSeg * seg,Double_t& x1, Double_t& y1, Double_t& z1, Double_t& x2, Double_t& y2, Double_t& z2)
{
    // calculates 2 coordinate points from segment

    Double_t teta = seg->getTheta();
    Double_t phi  = seg->getPhi();
    Double_t z    = seg->getZ();
    Double_t r    = seg->getR();
    Double_t pi2  = acos(-1.)/2.;

    Double_t X = r * cos(phi + pi2);
    Double_t Y = r * sin(phi + pi2);
    Double_t Z = z;

    x1 = X;
    y1 = Y;
    z1 = Z;
    x2 = X + cos(phi) * sin(teta);
    y2 = Y + sin(phi) * sin(teta);
    z2 = Z + cos(teta);

}

void  HMdcDeDx2::findBin(Int_t m,Double_t* angle,Double_t* mindist,Int_t* abin,Int_t* dbin)
{
    Double_t a =* angle;
    Double_t d =* mindist;

    if(d>MaxMdcMinDist[m]) d = MaxMdcMinDist[m];
    if(a > 90) a = 90.;
    if(a < 0)  a =  0.;

    (*abin) = (Int_t)((90. - a) / AngleStep);
    (*dbin) = (Int_t)(d / MinDistStep[m]);

    if((*abin) == 18)     (*abin) = 17;
    if((*dbin) == N_DIST) (*dbin) = N_DIST - 1;

}

Double_t HMdcDeDx2::toTSigma(Int_t s,Int_t m,Double_t angle,Double_t mindist,Int_t shift)
{
    // returns asymmetric width of ToT for single drift cell
    // shift == 0 (default) : returns 0
    //       ==-1           : low   bound width of ToT distribution (sigma)
    //       == 1           : upper bound width of ToT distribution (sigma)

    Int_t abin = 0;
    Int_t dbin = 0;

    findBin(m,&angle,&mindist,&abin,&dbin);

    // shift is half FWHM = 2.355/2=1.1775 sigma
    // we have to recalulate sigma here

    if      (shift < 0) return (shiftpar[s][m][abin][dbin][0] / 1.1775);
    else if (shift > 0) return (shiftpar[s][m][abin][dbin][1] / 1.1775);
    else return 0.;
}

Double_t HMdcDeDx2::toTTodEdX(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t ToT)
{
    // calculation of ToT -> dEdX for single drift cell

    Int_t abin = 0;
    Int_t dbin = 0;

    findBin(m,&angle,&mindist,&abin,&dbin);

    Double_t* p= &par[s][m][abin][dbin][0];
    //-----------------------------------------------------

    if(ToT > parMax[s][m][abin][dbin]) ToT = parMax[s][m][abin][dbin];

    Double_t dEdX = TMath::Power(10.,(TMath::Power((ToT - p[0]) / p[1],(1. / p[2])))) - p[3];
    if(!TMath::Finite(dEdX)){
	Warning("totTodEdX()","dEdX out of range: s %i m %i abin %i dbin %i ToT %1.15e dEdX %1.15e !",s,m,abin,dbin,ToT,dEdX);
	dEdX = -1;
    }
    return dEdX;
}
Double_t HMdcDeDx2::dEdXToToT(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t dEdX)
{
    // calculation of dEdX -> ToT for single drift cell

    Int_t abin = 0;
    Int_t dbin = 0;

    findBin(m,&angle,&mindist,&abin,&dbin);


    Double_t* p = &par[s][m][abin][dbin][0];

    Double_t ToT = p[0] + p[1] * TMath::Power((TMath::Log10(dEdX +  p[3])),p[2]);

    if(!TMath::Finite(ToT)){
	Warning("dEdXToToT()","ToT out of range: s %i m %i abin %i dbin %i ToT %1.15e dEdX %1.15e !",s,m,abin,dbin,ToT,dEdX);
	ToT = -1;
    }

    return ToT;
}
Double_t HMdcDeDx2::normalize(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t ToT)
{
    // calibrate ToT :
    // dEdX=TotTodEdX(t2-t1) ----> dEdXToToT(dEdX) for reference module,angle,distbin
    if(ToT < 0)  return ToT;
    Double_t   dEdX = toTTodEdX(s,m,angle,mindist,ToT);
    if(dEdX >= 0) ToT = dEdXToToT(s,ref_MOD,ref_ANGLE,ref_DIST,dEdX);

    return ToT;
}



//----------------------dEdx-----------------------------------------------

Double_t HMdcDeDx2::energyLoss(Int_t id,Double_t p,Double_t hefr)
{
    // Calculates the dEdx (MeV 1/g cm2) of an particle with GEANT ID id
    // and momentum p (MeV) for He/i-butan gas mixture with He fraction hefr
    // (He (hefr) / i-butan (1-hefr)). The fomular is taken from PDG and doesn't
    // include the density correction term. The values for the mean excitation
    // energy are taken from Sauli.

    if(p == 0)                 return -1;
    if(hefr < 0. || hefr > 1.) return -1;
    Double_t mass   = HPidPhysicsConstants::mass(id);
    if(mass == 0) return -1;

    // keep function value inside boundary
    if      (mass <  150                 && p < 5 )   p = 5 ;
    else if (mass >= 150  && mass < 1000 && p < 20)   p = 20;
    else if (mass >= 1000 && mass < 2000 && p < 35)   p = 35;
    else if (mass >= 2000                && p < 55)   p = 55;

    // Z and A
    Double_t Z_gas  = 2. * hefr + (1. - hefr) * 34.;
    Double_t A_gas  = 4. * hefr + (1. - hefr) * 58.;
    // I_0 and I
    Double_t I_0_gas = 24.6 * hefr + (1. - hefr) * 10.8;
    Double_t I2      = pow(I_0_gas * Z_gas * (1.e-6),2); // sauli
    //Double_t I2     =pow(16. * pow(Z_gas,0.9),2); //C.Lippmann thesis


    Double_t K      = 0.307075; // MeV cm2 PDG, 4 * pi * N_{A} * r_{e}2 * m_{e}2 * c2
    Double_t mass2  = pow(mass,2);
    Double_t m_ec2  = HPidPhysicsConstants::mass(3);
    Double_t z2     = pow((Float_t)HPidPhysicsConstants::charge(id),2);
    Double_t p2     = pow(p,2);
    Double_t beta2  = 1. / ((mass2/p2) + 1.);
    Double_t gamma2 = 1./ (1 - beta2);
    Double_t gamma  = sqrt(gamma2);

    Double_t Tmax   = (2. * m_ec2 * beta2 * gamma2) / (1. + 2.* gamma * (m_ec2 / mass) + pow((m_ec2 / mass),2));
    Double_t term1  = K * z2 * (Z_gas / A_gas) * (1. / beta2);
    Double_t term2  = ((2. * m_ec2 * beta2 * gamma2 * Tmax) / I2);
    Double_t dedx   = term1 * (0.5 * log(term2) - beta2);

    return dedx;
}
TGraph* HMdcDeDx2::energyLossGraph(Int_t id,Double_t hefr,TString opt,Bool_t exchange,Int_t markerstyle,Int_t markercolor,Float_t markersize)
{
    // creates a TGraph for particle with GEANT ID id and
    // and He/i-butan gas mixture with He fraction hefr
    // (He (hefr) / i-butan (1-hefr) dEdx vs p,beta,1/beta2 or betagamma
    // depending on opt (p,beta,1/beta2,betagamma). exchange=kTRUE
    // will exchange x-axis and y-axis.


    Double_t pmin = 50;
    Double_t pmax = 100000.;
    Int_t    np   = 100;

    Double_t ptot = 0;
    Double_t xarray[np];
    Double_t yarray[np];

    Int_t    vers = 0;
    Double_t mass = HPidPhysicsConstants::mass(id);


    if(opt.CompareTo("p"             ) == 0) vers = 0;
    else if(opt.CompareTo("beta"     ) == 0) vers = 1;
    else if(opt.CompareTo("1/beta2"  ) == 0) vers = 2;
    else if(opt.CompareTo("betagamma") == 0) vers = 3;
    else {cout<<"HMdcDedx::calcDeDxGraph():unknow option!"<<endl;}

    for(Int_t i = 1; i <= np; i ++)
    {
      ptot = pmin * pow((pmax / pmin),(i - 1) / (Double_t)(np - 1.));
      if(vers == 0){xarray[i - 1] = ptot;}
      if(vers == 1){xarray[i - 1] = sqrt(1. / ((pow(mass,2) / pow(ptot,2)) + 1.));}
      if(vers == 2){xarray[i-1] = ((pow(mass,2) / pow(ptot,2)) + 1.);}
      if(vers == 3){xarray[i-1] = (ptot / mass);}
      yarray[i - 1] = HMdcDeDx2::energyLoss(id,ptot,hefr);
    }
    Char_t name[300];
    sprintf(name,"dedx_%s_He_%5.1f_ibutane_%5.1f",HPidPhysicsConstants::pid(id),hefr * 100,(1 - hefr) * 100);
    TGraph* g = 0;
    if(!exchange)g = new TGraph(np,xarray,yarray);
    else         g = new TGraph(np,yarray,xarray);
    g->SetName(name);
    g->SetMarkerStyle(markerstyle);
    g->SetMarkerSize (markersize);
    g->SetMarkerColor(markercolor);

    return g;
}
Double_t HMdcDeDx2::scaledTimeAboveThreshold(HGeantKine* kine,
					     Double_t p,
                                             Float_t t1, Float_t t2, Float_t* t2err,
					     Int_t s,Int_t m,
					     Double_t alpha,Double_t mindist)
{
    // calculated the t2 of drift cell measurent
    // from dedx of a given particle with momentum p.
    // The assumption is that all drift cell measurements
    // belonging to one track can be treated independent.
    // return t2 and smeared error to pointer t2err for an
    // single measurement. If the result t2-t1 would
    // be < 0 a random value for t2 0-20 ns larger than
    // t1(already including error) will be created.
    // The smeared error will be returned to the pointer t2err.


    //------------------------------------------------
    // check input
    if(!kine) {
        Error("caledTimeAboveThreshold()","retrieved ZERO pointer for kine object!");
	return -99;
    }
    if(p < 0) {
        Error("caledTimeAboveThreshold()","momentum = %5.3f is wrong!",p);
	return -99;
    }
    if(t1 < -997) {
        Error("caledTimeAboveThreshold()","t1 = %5.3f is wrong!",t1);
	return -99;
    }
    if(t2 < -997) {
        Error("caledTimeAboveThreshold()","t2 = %5.3f is wrong!",t1);
	return -99;
    }
    //------------------------------------------------

    //------------------------------------------------
    // get particle ID from kine object
    Int_t pTr = -1,pID = 0;
    kine->getParticle(pTr,pID);
    //------------------------------------------------


    //------------------------------------------------
    // get mass and charge of particle for dedx
    // calculation
    Double_t mass   = HPidPhysicsConstants::mass  (pID);
    Int_t    charge = HPidPhysicsConstants::charge(pID);


    if(charge == 0 || mass == 0)
    {
        Warning("HMdcDeDx2::scaledEnergyLoss()","id %i %s mass %7.5f charge %i p %7.5f skipped!"
                ,pID,HPidPhysicsConstants::pid(pID),mass,charge,p);
        return t2;
    }
    //------------------------------------------------

    //------------------------------------------------
    // calulate analytical dedx of particle

    Double_t dedx   = energyLoss(pID,p,getHeFraction());

    if(!TMath::Finite(dedx)){
	Error("caledTimeAboveThreshold()","dedx either NaN or Inf!");
        return t2;
    }
    //------------------------------------------------

    //------------------------------------------------
    // recalculate dedx to mean value of Time over
    // threshold for the given drift chamber
    Double_t ToTmean = dEdXToToT(s,m,alpha,mindist,dedx);

    if(!TMath::Finite(ToTmean) || dedx < 0){
	Error("caledTimeAboveThreshold()","ToT either NaN or Inf!");
	cout<<"s "      <<s
            <<" m "     <<m
	    <<" angle " <<alpha
	    <<" dist "  <<mindist
	    <<" dedx "  <<dedx
	    <<" pid "   <<pID
	    <<" mass "  <<mass
	    <<" charge "<<charge
	    <<" mom "   <<p
	    <<endl;

	return t2;
    }
    //------------------------------------------------

    //------------------------------------------------
    // generate some realistic smearing arround
    // the ToT mean value for single measurement
    // The distribution of ToT is asymmetric.
    // the smearing has to be done for left and right half
    // different.
    Double_t sigL   = toTSigma (s,m,alpha,mindist,-1);
    Double_t sigR   = toTSigma (s,m,alpha,mindist, 1);

    Double_t smear;
    if(gRandom->Rndm() > 0.5){
	smear =  TMath::Abs(gRandom->Gaus(0., sigR));
    } else{
	smear = -TMath::Abs(gRandom->Gaus(0., sigL));
    }
    if(t2err) * t2err = smear;
    //------------------------------------------------


    //------------------------------------------------
    // return t2 mean of an  single measurement. If the
    // result t2-t1 would  be < 0 a random value for
    // t2 0-20 ns larger than t1 will be created.
    t2 = ToTmean + t1;
    if(t2 <= t1) t2 = t1 + gRandom->Rndm() * 20.;

    if(debug){
	cout<<"scaledEnergyLoss() "
	    <<" s "     <<s
            <<" m "     <<m
	    <<" a "     <<alpha
            <<" d "     <<mindist
	    <<" id "    <<pID
	    <<" mass "  <<mass
	    <<" charge "<<charge
	    <<" mom "   <<p
            <<" dedx "  <<dedx
	    <<" t1 "    <<t1
            <<" t2 "    <<t2
	    <<" t2err " <<smear
            <<" tot "   <<ToTmean
	    <<" sigL "  <<sigL
	    <<" sigR "  <<sigR
            <<endl;
    }

    return t2;
}