]>
Commit | Line | Data |
---|---|---|
284050f7 | 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 | } |