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