#include "TMath.h"
#include "TRandom.h"

#include <Riostream.h>

//______________________________________________________________________________
Double_t XMedian(Int_t n, const Double_t *arr)
{
   if (n == 1) return arr[0];

// Create index and sort array
//x   Int_t *index = 0;
//x   if (!(index = new Int_t[n])) {return 1;} //return not ok
   Long64_t *index = 0;
   if (!(index = new Long64_t[n])) {return 1;} //return not ok
   TMath::Sort(n, arr, index);

// Find median
   Int_t k;
   Double_t median = 0;
   if ((n % 2) == 0){
      k = (Int_t)TMath::Floor(n / 2.0) - 1;
      median = (arr[index[k]] + arr[index[k+1]])/2;
   } else {
      k = (Int_t)TMath::Floor(n / 2.0);
      median = arr[index[k]];
   }//if

   delete [] index;

   return median;
}//Median

//______________________________________________________________________________
Double_t *Weights(Int_t n, const Double_t *arr, Double_t *w)
{
   // Weighting function using Tukey biweight.

   Double_t *x   = new Double_t[n];
   Double_t  med = XMedian(n, arr);
   Double_t  max = 0.0;
   Double_t  sum = 0.0;

   for (Int_t i=0; i<n; i++) {
      x[i] = arr[i] - med;
      max  = (TMath::Abs(x[i]) > max) ? TMath::Abs(x[i]) : max;
   }//for_i

   if (max == 0.0) {
      for (Int_t i=0; i<n; i++) {
         w[i] = arr[i];
         sum += w[i];
      }//for_i
   } else {
      // Tukey biweight
      Double_t uu  = 0.0;
      Double_t uu2 = 0.0;
      for (Int_t i=0; i<n; i++) {
         uu   = x[i]/max;
         uu2  = (1.0 - uu*uu);
         w[i] = uu2*uu2;
         sum += w[i];
      }//for_i
   }//if

   // final weight
   sum = (sum == 0.0) ? n : sum;  //prevent zero division
   for (Int_t i=0; i<n; i++) {
      w[i] = w[i]/sum;
   }//for_i

   delete [] x;

   return w;
}//Weights


//______________________________________________________________________________
void TestWeights()
{
   Int_t n = 24;

   Double_t data[] = {2.665500, 5.037562, 4.501229, 5.010825, 
                      4.993828, 5.481654, 8.313473, 4.997974, 
                      2.949897, 2.586380, 3.554085, 7.800919, 
                      2.275865, 3.969712, 3.762094, 3.905486, 
                      2.761068, 2.820203, 5.120746, 3.908563, 
                      3.791251, 3.422354, 4.651827, 3.649799};

   Double_t *w = new Double_t[n];

   w = Weights(n, data, w);

   cout << "w:"; for (Int_t i=0;i<n;i++) cout << " - " << w[i]; cout << endl; 

   delete [] w;
}//TestMedian