]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliMultiDimVector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliMultiDimVector.cxx
CommitLineData
ac7636a0 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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
7b45817b 16/* $Id$ */
ac7636a0 17
18///////////////////////////////////////////////////////////////////
19// //
20// Implementation of the class to store the number of signal //
21// and background events in bins of the cut values //
22// Origin: Elena Bruna (bruna@to.infn.it) //
23// Updated: Sergey Senyukov (senyukov@to.infn.it) //
3aa6ce53 24// Francesco Prino (prino@to.infn.it) //
25// Last Updated: Giacomo Ortona (ortona@to.infn.it) //
ac7636a0 26// //
27///////////////////////////////////////////////////////////////////
28
0adeb69c 29#include <fstream>
30#include <Riostream.h>
ac7636a0 31#include "TH2.h"
32#include "AliMultiDimVector.h"
33#include "AliLog.h"
ac7636a0 34#include "TString.h"
35
c64cb1f6 36using std::cout;
37using std::endl;
38
ac7636a0 39ClassImp(AliMultiDimVector)
40//___________________________________________________________________________
41AliMultiDimVector::AliMultiDimVector():TNamed("AliMultiDimVector","default"),
42fNVariables(0),
43fNPtBins(0),
44fVett(0),
45fNTotCells(0),
46fIsIntegrated(0)
47{
48 // default constructor
57291f02 49
50 for(Int_t i=0; i<fgkMaxNVariables; i++) {
51 fNCutSteps[i]=0;
52 fMinLimits[i]=0;
53 fMaxLimits[i]=0;
54 fGreaterThan[i]=kTRUE;
55 fAxisTitles[i]="";
56 }
57 for(Int_t i=0; i<fgkMaxNPtBins+1; i++) {
58 fPtLimits[i]=0;
59 }
ac7636a0 60}
61//___________________________________________________________________________
2350a1d9 62AliMultiDimVector::AliMultiDimVector(const char *name,const char *title, const Int_t nptbins, const Float_t* ptlimits, const Int_t npars, const Int_t *nofcells,const Float_t *loosecuts, const Float_t *tightcuts, const TString *axisTitles):TNamed(name,title),
ac7636a0 63fNVariables(npars),
64fNPtBins(nptbins),
65fVett(0),
66fNTotCells(0),
67fIsIntegrated(0){
68// standard constructor
57291f02 69
70 for(Int_t i=0; i<fgkMaxNVariables; i++) {
71 fNCutSteps[i]=0;
72 fMinLimits[i]=0;
73 fMaxLimits[i]=0;
74 fGreaterThan[i]=kTRUE;
75 fAxisTitles[i]="";
76 }
77 for(Int_t i=0; i<fgkMaxNPtBins+1; i++) {
78 fPtLimits[i]=0;
79 }
80
ac7636a0 81 ULong64_t ntot=1;
82 for(Int_t i=0;i<fNVariables;i++){
83 ntot*=nofcells[i];
84 fNCutSteps[i]=nofcells[i];
3aa6ce53 85 if(loosecuts[i] == tightcuts[i]){
86 if(loosecuts[i]!=0){
87 printf("AliMultiDimVector::AliMultiDimVector: WARNING! same tight/loose variable for variable number %d. AliMultiDimVector with run with the following values: loose: %f; tight: %f\n",i,tightcuts[i]-0.1*tightcuts[i],tightcuts[i]);
88 fMinLimits[i]=tightcuts[i]-0.1*tightcuts[i];
89 fMaxLimits[i]=tightcuts[i];
90 }else{
91 fMinLimits[i]=0;
92 fMaxLimits[i]=0.0001;
93 }
94 fGreaterThan[i]=kTRUE;
95 }
96 if(loosecuts[i] < tightcuts[i]){
ac7636a0 97 fMinLimits[i]=loosecuts[i];
98 fMaxLimits[i]=tightcuts[i];
99 fGreaterThan[i]=kTRUE;
100 }else{
101 fMinLimits[i]=tightcuts[i];
102 fMaxLimits[i]=loosecuts[i];
103 fGreaterThan[i]=kFALSE;
104 }
105 fAxisTitles[i]=axisTitles[i].Data();
106 }
bac542b0 107 fNTotCells=ntot*fNPtBins;
ac7636a0 108 fVett.Set(fNTotCells);
bac542b0 109 for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=ptlimits[ipt];
110 for(Int_t ipt=fNPtBins+1;ipt<fgkMaxNPtBins+1;ipt++) fPtLimits[ipt]=999.;
ac7636a0 111 for (ULong64_t j=0;j<fNTotCells;j++) fVett.AddAt(0,j);
112}
113//___________________________________________________________________________
114AliMultiDimVector::AliMultiDimVector(const AliMultiDimVector &mv):TNamed(mv.GetName(),mv.GetTitle()),
115fNVariables(mv.fNVariables),
116fNPtBins(mv.fNPtBins),
117fVett(0),
118fNTotCells(mv.fNTotCells),
119fIsIntegrated(mv.fIsIntegrated)
120{
121 // copy constructor
57291f02 122
123 for(Int_t i=0; i<fgkMaxNVariables; i++) {
124 fNCutSteps[i]=0;
125 fMinLimits[i]=0;
126 fMaxLimits[i]=0;
127 fGreaterThan[i]=kTRUE;
128 fAxisTitles[i]="";
129 }
130 for(Int_t i=0; i<fgkMaxNPtBins+1; i++) {
131 fPtLimits[i]=0;
132 }
133
134
135 for(Int_t i=0;i<fNVariables;i++){
136 fNCutSteps[i]=mv.GetNCutSteps(i);
137 fMinLimits[i]=mv.GetMinLimit(i);
138 fMaxLimits[i]=mv.GetMaxLimit(i);
139 fGreaterThan[i]=mv.GetGreaterThan(i);
140 fAxisTitles[i]=mv.GetAxisTitle(i);
141 }
142 fVett.Set(fNTotCells);
143
144 for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv.GetPtLimit(ipt);
145 for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
146}
147//___________________________________________________________________________
148AliMultiDimVector &AliMultiDimVector::operator=(const AliMultiDimVector &mv)
149{
150 // assignment operator
151
152 if(&mv == this) return *this;
153
154 TNamed::operator=(mv);
155
156 fNVariables=mv.fNVariables;
157 fNPtBins=mv.fNPtBins;
158 fNTotCells=mv.fNTotCells;
159 fIsIntegrated=mv.fIsIntegrated;
160
161
162 for(Int_t i=0; i<fgkMaxNVariables; i++) {
163 fNCutSteps[i]=0;
164 fMinLimits[i]=0;
165 fMaxLimits[i]=0;
166 fGreaterThan[i]=kTRUE;
167 fAxisTitles[i]="";
168 }
169 for(Int_t i=0; i<fgkMaxNPtBins+1; i++) {
170 fPtLimits[i]=0;
171 }
172
173
ac7636a0 174 for(Int_t i=0;i<fNVariables;i++){
175 fNCutSteps[i]=mv.GetNCutSteps(i);
176 fMinLimits[i]=mv.GetMinLimit(i);
177 fMaxLimits[i]=mv.GetMaxLimit(i);
178 fGreaterThan[i]=mv.GetGreaterThan(i);
179 fAxisTitles[i]=mv.GetAxisTitle(i);
180 }
181 fVett.Set(fNTotCells);
bac542b0 182
183 for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv.GetPtLimit(ipt);
ac7636a0 184 for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
c5d30c64 185
186 return *this;
ac7636a0 187}
188//___________________________________________________________________________
189void AliMultiDimVector::CopyStructure(const AliMultiDimVector* mv){
190// Sets dimensions and limit from mv
191 fNVariables=mv->GetNVariables();
192 fNPtBins=mv->GetNPtBins();
193 fNTotCells=mv->GetNTotCells();
194 fIsIntegrated=mv->IsIntegrated();
195 for(Int_t i=0;i<fNVariables;i++){
196 fNCutSteps[i]=mv->GetNCutSteps(i);
197 fMinLimits[i]=mv->GetMinLimit(i);
198 fMaxLimits[i]=mv->GetMaxLimit(i);
199 fGreaterThan[i]=mv->GetGreaterThan(i);
200 fAxisTitles[i]=mv->GetAxisTitle(i);
201 }
bac542b0 202 for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv->GetPtLimit(ipt);
203 fVett.Set(fNTotCells);
ac7636a0 204}
205//______________________________________________________________________
206Bool_t AliMultiDimVector::GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const {
207// returns matrix element indices and Pt bin from global index
208 if(globadd>=fNTotCells) return kFALSE;
209 ULong64_t r=globadd;
210 Int_t prod=1;
211 Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
212 for(Int_t k=0;k<fNVariables;k++) nOfCellsPlusLevel[k]=fNCutSteps[k];
213 nOfCellsPlusLevel[fNVariables]=fNPtBins;
214
215 for(Int_t i=0;i<fNVariables+1;i++) prod*=nOfCellsPlusLevel[i];
216 for(Int_t i=0;i<fNVariables+1;i++){
217 prod/=nOfCellsPlusLevel[i];
218 if(i<fNVariables) ind[i]=r/prod;
219 else ptbin=r/prod;
220 r=globadd%prod;
221 }
222 return kTRUE;
223}
224//______________________________________________________________________
225Bool_t AliMultiDimVector::GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const {
226 Int_t ind[fgkMaxNVariables];
227 Bool_t retcode=GetIndicesFromGlobalAddress(globadd,ind,ptbin);
228 if(!retcode) return kFALSE;
229 for(Int_t i=0;i<fNVariables;i++) cuts[i]=GetCutValue(i,ind[i]);
230 return kTRUE;
231}
232//______________________________________________________________________
2350a1d9 233ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(const Int_t *ind, Int_t ptbin) const {
ac7636a0 234 // Returns the global index of the cell in the matrix
235 Int_t prod=1;
236 ULong64_t elem=0;
237 Int_t indexPlusLevel[fgkMaxNVariables+1];
238 Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
239 for(Int_t i=0;i<fNVariables;i++){
240 indexPlusLevel[i]=ind[i];
241 nOfCellsPlusLevel[i]=fNCutSteps[i];
242 }
243 indexPlusLevel[fNVariables]=ptbin;
244 nOfCellsPlusLevel[fNVariables]=fNPtBins;
245
246 for(Int_t i=0;i<fNVariables+1;i++){
247 prod=indexPlusLevel[i];
248 if(i<fNVariables){
249 for(Int_t j=i+1;j<fNVariables+1;j++){
250 prod*=nOfCellsPlusLevel[j];
251 }
252 }
253 elem+=prod;
254 }
255 return elem;
256}
257//______________________________________________________________________
2350a1d9 258Bool_t AliMultiDimVector::GetIndicesFromValues(const Float_t *values, Int_t *ind) const {
ac7636a0 259 // Fills the array of matrix indices strating from variable values
260 for(Int_t i=0;i<fNVariables;i++){
261 if(fGreaterThan[i]){
262 if(values[i]<GetMinLimit(i)) return kFALSE;
263 ind[i]=(Int_t)((values[i]-fMinLimits[i])/GetCutStep(i));
264 if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
265 }else{
266 if(values[i]>GetMaxLimit(i)) return kFALSE;
267 ind[i]=(Int_t)((fMaxLimits[i]-values[i])/GetCutStep(i));
268 if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
269 }
270 }
271 return kTRUE;
272}
273//______________________________________________________________________
2350a1d9 274ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(const Float_t *values, Int_t ptbin) const {
ac7636a0 275 // Returns the global index of the cell in the matrix
276 Int_t ind[fgkMaxNVariables];
277 Bool_t retcode=GetIndicesFromValues(values,ind);
278 if(retcode) return GetGlobalAddressFromIndices(ind,ptbin);
279 else{
280 AliError("Values out of range");
281 return fNTotCells+999;
282 }
283}
284//_____________________________________________________________________________
285void AliMultiDimVector::MultiplyBy(Float_t factor){
286 // multiply the AliMultiDimVector by a constant factor
287 for(ULong64_t i=0;i<fNTotCells;i++){
2350a1d9 288 if(fVett.At(i)>0.)
ac7636a0 289 fVett.AddAt(fVett.At(i)*factor,i);
290 else fVett.AddAt(-1,i);
291 }
292
293}
294//_____________________________________________________________________________
2350a1d9 295void AliMultiDimVector::Multiply(const AliMultiDimVector* mv,Float_t factor){
ac7636a0 296 // Sets AliMultiDimVector=mv*constant factor
297 for(ULong64_t i=0;i<fNTotCells;i++){
2350a1d9 298 if(mv->GetElement(i)>0.)
ac7636a0 299 fVett.AddAt(mv->GetElement(i)*factor,i);
300 else fVett.AddAt(-1,i);
301 }
302}
303//_____________________________________________________________________________
2350a1d9 304void AliMultiDimVector::Multiply(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
ac7636a0 305 // Sets AliMultiDimVector=mv1*mv2
306 for(ULong64_t i=0;i<fNTotCells;i++){
2350a1d9 307 if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
ac7636a0 308 fVett.AddAt(mv1->GetElement(i)*mv2->GetElement(i),i);
309 else fVett.AddAt(-1,i);
310 }
311}
312//_____________________________________________________________________________
2350a1d9 313void AliMultiDimVector::Add(const AliMultiDimVector* mv){
ac7636a0 314 // Sums contents of mv to AliMultiDimVector
315 if (mv->GetNTotCells()!=fNTotCells){
316 AliError("Different dimension of the vectors!!");
317 }else{
318 for(ULong64_t i=0;i<fNTotCells;i++)
2350a1d9 319 if(mv->GetElement(i)>0. && fVett.At(i)>0.)
ac7636a0 320 fVett.AddAt(fVett.At(i)+mv->GetElement(i),i);
321 else fVett.AddAt(-1,i);
322 }
323}
324//_____________________________________________________________________________
2350a1d9 325void AliMultiDimVector::Sum(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
ac7636a0 326 // Sets AliMultiDimVector=mv1+mv2
327 if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
328 AliError("Different dimension of the vectors!!");
329 }
330 else{
331 for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
2350a1d9 332 if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
ac7636a0 333 fVett.AddAt(mv1->GetElement(i)+mv2->GetElement(i),i);
334 else fVett.AddAt(-1,i);
335 }
336 }
337}
338//_____________________________________________________________________________
2350a1d9 339void AliMultiDimVector::LinearComb(const AliMultiDimVector* mv1, Float_t norm1, const AliMultiDimVector* mv2, Float_t norm2){
ac7636a0 340 // Sets AliMultiDimVector=n1*mv1+n2*mv2
341 if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
342 AliError("Different dimension of the vectors!!");
343 }
344 else{
345 for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
2350a1d9 346 if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
ac7636a0 347 fVett.AddAt(norm1*mv1->GetElement(i)+norm2*mv2->GetElement(i),i);
348 else fVett.AddAt(-1,i);
349 }
350 }
351}
352//_____________________________________________________________________________
2350a1d9 353void AliMultiDimVector::DivideBy(const AliMultiDimVector* mv){
ac7636a0 354 // Divide AliMulivector by mv
355 if (mv->GetNTotCells()!=fNTotCells) {
356 AliError("Different dimension of the vectors!!");
357 }
358 else{
359 for(ULong64_t i=0;i<fNTotCells;i++)
2350a1d9 360 if(mv->GetElement(i)!=0 &&mv->GetElement(i)>0. && fVett.At(i)>0.)
ac7636a0 361 fVett.AddAt(fVett.At(i)/mv->GetElement(i),i);
362 else fVett.AddAt(-1,i);
363 }
364
365}
366//_____________________________________________________________________________
2350a1d9 367void AliMultiDimVector::Divide(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
ac7636a0 368 // Sets AliMultiDimVector=mv1/mv2
369 if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
370 AliError("Different dimension of the vectors!!");
371 }
372 else{
373 for(ULong64_t i=0;i<mv1->GetNTotCells();i++)
2350a1d9 374 if(mv2->GetElement(i)!=0&& mv2->GetElement(i)>0.&& mv1->GetElement(i)>0.)
ac7636a0 375 {
376 fVett.AddAt(mv1->GetElement(i)/mv2->GetElement(i),i);
377 }
378 else fVett.AddAt(-1,i);
379 }
380}
381//_____________________________________________________________________________
382void AliMultiDimVector::Sqrt(){
383 // Sqrt of elements of AliMultiDimVector
384 for(ULong64_t i=0;i<fNTotCells;i++) {
385 if(fVett.At(i)>=0) fVett.AddAt(TMath::Sqrt(fVett.At(i)),i);
386 else {
387 fVett.AddAt(-1,i);
388 }
389 }
390}
391//_____________________________________________________________________________
2350a1d9 392void AliMultiDimVector::Sqrt(const AliMultiDimVector* mv){
ac7636a0 393 // Sets AliMultiDimVector=sqrt(mv)
394 for(ULong64_t i=0;i<fNTotCells;i++)
395 if(mv->GetElement(i)>=0) fVett.AddAt(TMath::Sqrt(mv->GetElement(i)),i);
396 else fVett.AddAt(-1,i);
397}
398//_____________________________________________________________________________
399void AliMultiDimVector::FindMaximum(Float_t& maxValue, Int_t *ind , Int_t ptbin){
400 // finds the element with maximum contents
401 const ULong64_t nelem=fNTotCells/fNPtBins;
402 TArrayF vett;
403 vett.Set(nelem);
404 ULong64_t runningAddress;
405 for(ULong64_t i=0;i<nelem;i++){
406 runningAddress=ptbin+i*fNPtBins;
407 vett.AddAt(fVett[runningAddress],i);
408 }
409 maxValue=TMath::MaxElement(nelem,vett.GetArray());
410 ULong64_t maxAddress=TMath::LocMax(nelem,vett.GetArray());
411 ULong64_t maxGlobalAddress=ptbin+maxAddress*fNPtBins;
412 Int_t checkedptbin;
413 GetIndicesFromGlobalAddress(maxGlobalAddress,ind,checkedptbin);
414}
415
0adeb69c 416//_____________________________________________________________________________
417//Int_t* AliMultiDimVector::FindLocalMaximum(Float_t& maxValue, Bool_t *isFree,Int_t* indFixed, Int_t ptbin){
418Int_t* AliMultiDimVector::FindLocalMaximum(Float_t& maxValue, Int_t *numFixed,Int_t* indFixed, Int_t nfixed,Int_t ptbin){
419 //return the elements with maximum content (maxValue) given fixed step for not free variables
420 //numFixed[nfixed] is the indices of the fixed variables in the cuts array [fNVariables]={kTRUE,kTRUE,...,kFALSE,...,kFALSE,...,kTRUE}
421 //indFixed[nfixed]={1,2} //nfixed is the number of false in isFree; indFixed contains the step for the i-th variable
422 //!!take care of deleting the array of index returned!!
423
424 // Int_t nfixed=0,nfree=0;
425 //Int_t indtmp[fNVariables];
426 if(nfixed>fNVariables)cout<<"AliMultiDimVector::FindLocalMaximum:ERROR! too many variables"<<endl;
427 ULong64_t nelem=1;
428 Int_t* indMax=new Int_t[fNVariables];
429 //Get the number of fixed vars
430 /*
431 for (Int_t iv=0;iv<fNVariables;iv++){
432 if(isFree[iv]){
433 nfree++;
434 nelem*=fNCutSteps[iv];
435 indMax[iv]=0;
436 }
437 else {
438 indMax[iv]=indFixed[nfixed];
439 if(indFixed[nfixed]>=GetNCutSteps(iv)){
440 indMax[iv]=0;
441 cout<<"AliMultiDimVector::FindLocalMaximum:ERROR! called fixed ind "<< indFixed[nfixed]<<" but "<<iv<<" var has only "<<GetNCutSteps(iv)<<" steps"<<endl;
442 }
443 nfixed++;
444 }
445 }
446 */
447 for (Int_t iv=0;iv<fNVariables;iv++)indMax[iv]=0;
448 for(Int_t i=0;i<nfixed;i++){
449 indMax[numFixed[i]]=indFixed[i];
450 }
451 //Get position of fixed vars
452 /*
453 Int_t fixedIndexes[nfixed];
454 Int_t iforfixed=0;
455 for (Int_t iv=0;iv<fNVariables;iv++){
456 if(!isFree[iv]){
457 fixedIndexes[iforfixed]=iv;
458 iforfixed++;
459 }
460 }
461 */
462 TArrayF vett;
463 vett.Set(nelem);
464
465 ULong64_t first=fNTotCells/fNPtBins*ptbin;
466 ULong64_t last=first+fNTotCells/fNPtBins;
467 Int_t dummyptbin;
468
469 maxValue=fVett[GetGlobalAddressFromIndices(indMax,ptbin)];
470 Int_t tmpInd[fNVariables];
471
472 //loop on multidimvector global addresses
473 for(ULong64_t iga=first;iga<last;iga++){
474 GetIndicesFromGlobalAddress(iga,tmpInd,dummyptbin);
475 Bool_t goodCell=kTRUE;
476 for(Int_t ifix=0;ifix<nfixed&&goodCell;ifix++){
477 // if(indFixed[ifix]!=indMax[fixedIndexes[ifix]])goodCell=kFALSE;
478 if(indFixed[ifix]!=tmpInd[numFixed[ifix]])goodCell=kFALSE;
479 }
480 if(goodCell){
481 if(fVett[iga]>maxValue){
482 maxValue=fVett[iga];
483 // GetIndicesFromGlobalAddress(iga,indMax,dummyptbin);
484 for(Int_t inv=0;inv<fNVariables;inv++)indMax[inv]=tmpInd[inv];
485 }
486 }
487 }
488
489 return indMax;
490}
491
ac7636a0 492//_____________________________________________________________________________
2350a1d9 493TH2F* AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, const Int_t* fixedVars, Int_t ptbin, Float_t norm){
ac7636a0 494 // Project the AliMultiDimVector on a 2D histogram
495
496 TString hisName=Form("hproj%s%dv%d",GetName(),secondVar,firstVar);
497 TString hisTit=Form("%s vs. %s",fAxisTitles[secondVar].Data(),fAxisTitles[firstVar].Data());
498 TH2F* h2=new TH2F(hisName.Data(),hisTit.Data(),fNCutSteps[firstVar],fMinLimits[firstVar],fMaxLimits[firstVar],fNCutSteps[secondVar],fMinLimits[secondVar],fMaxLimits[secondVar]);
499
500 Int_t index[fgkMaxNVariables];
501 for(Int_t i=0;i<fNVariables;i++){
502 index[i]=fixedVars[i];
503 }
504
505 for(Int_t i=0;i<fNCutSteps[firstVar];i++){
506 for(Int_t j=0;j<fNCutSteps[secondVar];j++){
507 index[firstVar]=i;
508 index[secondVar]=j;
509 Float_t cont=GetElement(index,ptbin)/norm;
510 Int_t bin1=i+1;
511 if(!fGreaterThan[firstVar]) bin1=fNCutSteps[firstVar]-i;
512 Int_t bin2=j+1;
513 if(!fGreaterThan[secondVar]) bin2=fNCutSteps[secondVar]-j;
514 h2->SetBinContent(bin1,bin2,cont);
515 }
516 }
517 return h2;
518}
519//_____________________________________________________________________________
520void AliMultiDimVector::GetIntegrationLimits(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
521 // computes bin limits for integrating the AliMultiDimVector
522 minbin=0;
523 maxbin=0;
524 if(iVar<fNVariables){
525 minbin=iCell;
526 maxbin=fNCutSteps[iVar]-1;
527 }
528}
529//_____________________________________________________________________________
530void AliMultiDimVector::GetFillRange(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
531 // computes range of cells passing the cuts for FillAndIntegrate
532 minbin=0;
533 maxbin=0;
534 if(iVar<fNVariables){
535 minbin=0; // bin 0 corresponds to loose cuts
536 maxbin=iCell;
537 }
538}
539//_____________________________________________________________________________
540void AliMultiDimVector::Integrate(){
541 // integrates the matrix
542 if(fIsIntegrated){
543 AliError("MultiDimVector already integrated");
544 return;
545 }
546 TArrayF integral(fNTotCells);
547 for(ULong64_t i=0;i<fNTotCells;i++) integral[i]=CountsAboveCell(i);
548 for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]= integral[i];
549 fIsIntegrated=kTRUE;
bac542b0 550}//_____________________________________________________________________________
2350a1d9 551ULong64_t* AliMultiDimVector::GetGlobalAddressesAboveCuts(const Float_t *values, Int_t ptbin, Int_t& nVals) const{
bac542b0 552 // fills an array with global addresses of cells passing the cuts
553
554 Int_t ind[fgkMaxNVariables];
555 Bool_t retcode=GetIndicesFromValues(values,ind);
556 if(!retcode){
557 nVals=0;
558 return 0x0;
559 }
560 for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
561 Int_t mink[fgkMaxNVariables];
562 Int_t maxk[fgkMaxNVariables];
563 Int_t size=1;
564 for(Int_t i=0;i<fgkMaxNVariables;i++){
565 GetFillRange(i,ind[i],mink[i],maxk[i]);
566 size*=(maxk[i]-mink[i]+1);
567 }
568 ULong64_t* indexes=new ULong64_t[size];
569 nVals=0;
570 for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
571 for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
572 for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
573 for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
574 for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
575 for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
576 for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
577 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
578 for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
579 for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
580 Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
581 indexes[nVals++]=GetGlobalAddressFromIndices(currentBin,ptbin);
582 }
583 }
584 }
585 }
586 }
587 }
588 }
589 }
590 }
591 }
592 return indexes;
ac7636a0 593}
594//_____________________________________________________________________________
595Float_t AliMultiDimVector::CountsAboveCell(ULong64_t globadd) const{
596 // integrates the counts of cells above cell with address globadd
597 Int_t ind[fgkMaxNVariables];
598 Int_t ptbin;
599 GetIndicesFromGlobalAddress(globadd,ind,ptbin);
600 for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
601 Int_t mink[fgkMaxNVariables];
602 Int_t maxk[fgkMaxNVariables];
603 for(Int_t i=0;i<fgkMaxNVariables;i++){
604 GetIntegrationLimits(i,ind[i],mink[i],maxk[i]);
605 }
606 Float_t sumcont=0.;
607 for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
608 for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
609 for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
610 for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
611 for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
612 for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
613 for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
614 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
615 for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
616 for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
617 Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
618 sumcont+=GetElement(currentBin,ptbin);
619 }
620 }
621 }
622 }
623 }
624 }
625 }
626 }
627 }
628 }
629 return sumcont;
630}
631//_____________________________________________________________________________
632void AliMultiDimVector::Fill(Float_t* values, Int_t ptbin){
633 // fills the cells of AliMultiDimVector corresponding to values
634 if(fIsIntegrated){
635 AliError("MultiDimVector already integrated -- Use FillAndIntegrate");
636 return;
637 }
638 Int_t ind[fgkMaxNVariables];
639 Bool_t retcode=GetIndicesFromValues(values,ind);
640 for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
641 if(retcode) IncrementElement(ind,ptbin);
642}
643//_____________________________________________________________________________
644void AliMultiDimVector::FillAndIntegrate(Float_t* values, Int_t ptbin){
645 // fills the cells of AliMultiDimVector passing the cuts
646 // The number of nested loops must match fgkMaxNVariables!!!!
647 fIsIntegrated=kTRUE;
648 Int_t ind[fgkMaxNVariables];
649 Bool_t retcode=GetIndicesFromValues(values,ind);
650 if(!retcode) return;
651 for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
652 Int_t mink[fgkMaxNVariables];
653 Int_t maxk[fgkMaxNVariables];
654 for(Int_t i=0;i<fgkMaxNVariables;i++){
655 GetFillRange(i,ind[i],mink[i],maxk[i]);
656 }
657 for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
658 for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
659 for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
660 for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
661 for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
662 for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
663 for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
664 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
665 for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
666 for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
667 Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
668 IncrementElement(currentBin,ptbin);
669 }
670 }
671 }
672 }
673 }
674 }
675 }
676 }
677 }
678 }
679
680}
681//_____________________________________________________________________________
2350a1d9 682void AliMultiDimVector::SuppressZeroBKGEffect(const AliMultiDimVector* mvBKG){
ac7636a0 683 // Sets to zero elements for which mvBKG=0
684 for(ULong64_t i=0;i<fNTotCells;i++)
2350a1d9 685 if(mvBKG->GetElement(i)<0.00000001) fVett.AddAt(0,i);
ac7636a0 686}
687//_____________________________________________________________________________
688AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBin){
689 // sums the elements of pt bins between firstBin and lastBin
690 if(firstBin<0 || lastBin>=fNPtBins || firstBin>=lastBin){
691 AliError("Bad numbers of Pt bins to be shrinked");
692 return 0;
693 }
694 Int_t nofcells[fgkMaxNVariables];
695 Float_t loosecuts[fgkMaxNVariables];
696 Float_t tightcuts[fgkMaxNVariables];
697 TString axisTitles[fgkMaxNVariables];
e11ae259 698 for(Int_t j=0;j<fgkMaxNVariables;j++) {
699 nofcells[j]=0;
700 loosecuts[j]=0.;
701 tightcuts[j]=0.;
702 axisTitles[j]="";
703 }
ac7636a0 704 for(Int_t i=0;i<fNVariables;i++){
705 nofcells[i]=fNCutSteps[i];
706 if(fGreaterThan[i]){
707 loosecuts[i]=fMinLimits[i];
708 tightcuts[i]=fMaxLimits[i];
709 }else{
710 loosecuts[i]=fMaxLimits[i];
711 tightcuts[i]=fMinLimits[i];
712 }
713 axisTitles[i]=fAxisTitles[i];
714 }
715 Int_t newNptbins=fNPtBins-(lastBin-firstBin);
bac542b0 716 Float_t ptlimits[fgkMaxNPtBins+1];
717 for(Int_t ipt=0; ipt<=firstBin;ipt++) ptlimits[ipt]=fPtLimits[ipt];
718 for(Int_t ipt=firstBin+1; ipt<newNptbins+1;ipt++) ptlimits[ipt]=fPtLimits[ipt+(lastBin-firstBin)];
719 AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),newNptbins,ptlimits,fNVariables,nofcells,loosecuts,tightcuts,axisTitles);
ac7636a0 720
721 ULong64_t nOfPointsPerPtbin=fNTotCells/fNPtBins;
722 ULong64_t addressOld,addressNew;
723 Int_t npb,opb;
724 for(npb=0;npb<firstBin;npb++){
725 opb=npb;
726 for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
727 addressOld=opb+k*fNPtBins;
728 addressNew=npb+k*newNptbins;
729 shrinkedMV->SetElement(addressNew,fVett[addressOld]);
730 }
731 }
732 npb=firstBin;
733 for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
734 Float_t summedValue=0.;
735 for(opb=firstBin;opb<=lastBin;opb++){
736 addressOld=opb+k*fNPtBins;
737 summedValue+=fVett[addressOld];
738 }
739 addressNew=npb+k*newNptbins;
740 shrinkedMV->SetElement(addressNew,summedValue);
741 }
742 for(npb=firstBin+1;npb<newNptbins;npb++){
743 opb=npb+(lastBin-firstBin);
744 for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
745 addressOld=opb+k*fNPtBins;
746 addressNew=npb+k*newNptbins;
747 shrinkedMV->SetElement(addressNew,fVett[addressOld]);
748 }
749 }
750 return shrinkedMV;
751}
bac542b0 752//_____________________________________________________________________________
3aa6ce53 753void AliMultiDimVector::SetNewLimits(Float_t* loose,Float_t* tight){
754 for(Int_t i=0;i<fNVariables;i++){
755 if(loose[i] < tight[i]){
756 fMinLimits[i]=loose[i];
757 fMaxLimits[i]=tight[i];
758 fGreaterThan[i]=kTRUE;
759 }else{
760 fMinLimits[i]=tight[i];
761 fMaxLimits[i]=loose[i];
762 fGreaterThan[i]=kFALSE;
763 }
764 }
765}
766//_____________________________________________________________________________
767void AliMultiDimVector::SwapLimits(Int_t ivar){
768 Float_t oldmin = fMinLimits[ivar];
769 fMinLimits[ivar] = fMaxLimits[ivar];
770 fMaxLimits[ivar] = oldmin;
771 if(fGreaterThan[ivar])fGreaterThan[ivar]=kFALSE;
772 else fGreaterThan[ivar]=kTRUE;
773}
774//_____________________________________________________________________________
bac542b0 775void AliMultiDimVector::PrintStatus(){
776 //
777 printf("Number of Pt bins = %d\n",fNPtBins);
778 printf("Limits of Pt bins = ");
779 for(Int_t ib=0;ib<fNPtBins+1;ib++) printf("%6.2f ",fPtLimits[ib]);
780 printf("\n");
781 printf("Number of cut variables = %d\n",fNVariables);
782 for(Int_t iv=0;iv<fNVariables;iv++){
783 printf("- Variable %d: %s\n",iv,fAxisTitles[iv].Data());
784 printf(" Nsteps= %d Rage = %6.2f %6.2f\n",
785 fNCutSteps[iv],fMinLimits[iv],fMaxLimits[iv]);
786 }
787}