STEER/AliMathBase.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 // Class AliMathBase
  19 //
  20 // Subset of  matheamtical functions  not included in the TMath
  21 //
  22
  23 ///////////////////////////////////////////////////////////////////////////
  24 #include "TMath.h"
  25 #include "AliMathBase.h"
  26 #include "Riostream.h"
  27
  28 ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
  29
  30 AliMathBase::AliMathBase() : TObject()
  31 {
  32 // Default constructor
  33 }
  34 ///////////////////////////////////////////////////////////////////////////
  35 AliMathBase::~AliMathBase()
  36 {
  37 // Destructor
  38 }
  39
  40
  41 //_____________________________________________________________________________
  42 void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
  43                            , Double_t &sigma, Int_t hh)
  44 {
  45   //
  46   // Robust estimator in 1D case MI version - (faster than ROOT version)
  47   //
  48   // For the univariate case
  49   // estimates of location and scatter are returned in mean and sigma parameters
  50   // the algorithm works on the same principle as in multivariate case -
  51   // it finds a subset of size hh with smallest sigma, and then returns mean and
  52   // sigma of this subset
  53   //
  54
  55   if (hh==0)
  56     hh=(nvectors+2)/2;
  57   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
  58   Int_t *index=new Int_t[nvectors];
  59   TMath::Sort(nvectors, data, index, kFALSE);
  60
  61   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  62   Double_t factor = faclts[nquant-1];
  63
  64   Double_t sumx  =0;
  65   Double_t sumx2 =0;
  66   Int_t    bestindex = -1;
  67   Double_t bestmean  = 0;
  68   Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
  69   for (Int_t i=0; i<hh; i++){
  70     sumx  += data[index[i]];
  71     sumx2 += data[index[i]]*data[index[i]];
  72   }
  73
  74   Double_t norm = 1./Double_t(hh);
  75   Double_t norm2 = 1./Double_t(hh-1);
  76   for (Int_t i=hh; i<nvectors; i++){
  77     Double_t cmean  = sumx*norm;
  78     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
  79     if (csigma<bestsigma){
  80       bestmean  = cmean;
  81       bestsigma = csigma;
  82       bestindex = i-hh;
  83     }
  84
  85     sumx  += data[index[i]]-data[index[i-hh]];
  86     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
  87   }
  88
  89   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
  90   mean  = bestmean;
  91   sigma = bstd;
  92   delete [] index;
  93
  94 }
  95
  96
  97
  98 void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh,  Float_t externalfactor)
  99 {
 100   // Modified version of ROOT robust EvaluateUni
 101   // robust estimator in 1D case MI version
 102   // added external factor to include precision of external measurement
 103   //
 104
 105   if (hh==0)
 106     hh=(nvectors+2)/2;
 107   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
 108   Int_t *index=new Int_t[nvectors];
 109   TMath::Sort(nvectors, data, index, kFALSE);
 110   //
 111   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
 112   Double_t factor = faclts[0];
 113   if (nquant>0){
 114     // fix proper normalization - Anja
 115     factor = faclts[nquant-1];
 116   }
 117
 118   //
 119   //
 120   Double_t sumx  =0;
 121   Double_t sumx2 =0;
 122   Int_t    bestindex = -1;
 123   Double_t bestmean  = 0;
 124   Double_t bestsigma = -1;
 125   for (Int_t i=0; i<hh; i++){
 126     sumx  += data[index[i]];
 127     sumx2 += data[index[i]]*data[index[i]];
 128   }
 129   //
 130   Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
 131   Double_t norm = 1./Double_t(hh);
 132   for (Int_t i=hh; i<nvectors; i++){
 133     Double_t cmean  = sumx*norm;
 134     Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
 135     if (csigma<bestsigma ||  bestsigma<0){
 136       bestmean  = cmean;
 137       bestsigma = csigma;
 138       bestindex = i-hh;
 139     }
 140     //
 141     //
 142     sumx  += data[index[i]]-data[index[i-hh]];
 143     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 144   }
 145
 146   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 147   mean  = bestmean;
 148   sigma = bstd;
 149   delete [] index;
 150 }
 151
 152
 153 //_____________________________________________________________________________
 154 Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
 155                         , Int_t *outlist, Bool_t down)
 156 {
 157   //
 158   //  Sort eleements according occurancy
 159   //  The size of output array has is 2*n
 160   //
 161
 162   Int_t * sindexS = new Int_t[n];     // temp array for sorting
 163   Int_t * sindexF = new Int_t[2*n];
 164   for (Int_t i=0;i<n;i++) sindexF[i]=0;
 165   //
 166   TMath::Sort(n,inlist, sindexS, down);
 167   Int_t last      = inlist[sindexS[0]];
 168   Int_t val       = last;
 169   sindexF[0]      = 1;
 170   sindexF[0+n]    = last;
 171   Int_t countPos  = 0;
 172   //
 173   //  find frequency
 174   for(Int_t i=1;i<n; i++){
 175     val = inlist[sindexS[i]];
 176     if (last == val)   sindexF[countPos]++;
 177     else{
 178       countPos++;
 179       sindexF[countPos+n] = val;
 180       sindexF[countPos]++;
 181       last =val;
 182     }
 183   }
 184   if (last==val) countPos++;
 185   // sort according frequency
 186   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
 187   for (Int_t i=0;i<countPos;i++){
 188     outlist[2*i  ] = sindexF[sindexS[i]+n];
 189     outlist[2*i+1] = sindexF[sindexS[i]];
 190   }
 191   delete [] sindexS;
 192   delete [] sindexF;
 193
 194   return countPos;
 195
 196 }