Temporary mod for CMake
[u/mrichter/AliRoot.git] / CORRFW / AliCFGridSparse.cxx
CommitLineData
db6722a5 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// AliCFGridSparse Class //
18// Class to accumulate data on an N-dimensional grid, to be used //
19// as input to get corrections for Reconstruction & Trigger efficiency//
20// Based on root THnSparse //
21// -- Author : S.Arcelli //
22// Still to be done: //
23// --Interpolate among bins in a range //
24//--------------------------------------------------------------------//
25//
26//
27#include "AliCFGridSparse.h"
28#include "THnSparse.h"
29#include "AliLog.h"
30#include "TMath.h"
31#include "TROOT.h"
32#include "TH1D.h"
33#include "TH2D.h"
34#include "TH3D.h"
35#include "TAxis.h"
36
37//____________________________________________________________________
38ClassImp(AliCFGridSparse)
39
40//____________________________________________________________________
41AliCFGridSparse::AliCFGridSparse() :
42 AliCFVGrid(),
db6722a5 43 fData(0x0)
44{
45 // default constructor
46}
47//____________________________________________________________________
48AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
49 AliCFVGrid(name,title),
db6722a5 50 fData(0x0)
51{
52 // default constructor
53}
54//____________________________________________________________________
55AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, const Int_t nVarIn, const Int_t * nBinIn, const Double_t *binLimitsIn) :
56 AliCFVGrid(name,title,nVarIn,nBinIn,binLimitsIn),
db6722a5 57 fData(0x0)
58{
59 //
60 // main constructor
61 //
62
63 fData=new THnSparseF(name,title,fNVar,fNVarBins);
64
65 if(binLimitsIn){
66 for(Int_t ivar=0;ivar<fNVar;ivar++){
67 Int_t nbins=fNVarBins[ivar]+1;
68 Double_t *array= new Double_t[nbins];
69 for(Int_t i=0;i<nbins;i++){
70 array[i]=fVarBinLimits[fOffset[ivar]+i];
71 }
72 fData->SetBinEdges(ivar, array);
73 delete array;
74 }
75 }
76}
77//____________________________________________________________________
78AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
79 AliCFVGrid(c),
db6722a5 80 fData(c.fData)
81{
82 //
83 // copy constructor
84 //
85 ((AliCFGridSparse &)c).Copy(*this);
86}
87
88//____________________________________________________________________
89AliCFGridSparse::~AliCFGridSparse()
90{
91 //
92 // destructor
93 //
94 if(fData) delete fData;
95}
96
97//____________________________________________________________________
98AliCFGridSparse &AliCFGridSparse::operator=(const AliCFGridSparse &c)
99{
100 //
101 // assigment operator
102 //
103 if (this != &c)
104 ((AliCFGridSparse &) c).Copy(*this);
105 return *this;
106}
107
108//____________________________________________________________________
109void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t *array)
110{
111 //
112 // setting the arrays containing the bin limits
113 //
114 fData->SetBinEdges(ivar, array);
115 //then fill the appropriate array in ALICFFrame, to be able to use
116 //the getter, in case....
117 Int_t nbins=fNVarBins[ivar]+1;
118 for(Int_t i=0;i<nbins;i++){
119 fVarBinLimits[fOffset[ivar]+i] =array[i];
120 }
121}
122
123//____________________________________________________________________
124void AliCFGridSparse::Fill(Double_t *var, Double_t weight)
125{
126 //
127 // Fill the grid,
128 // given a set of values of the input variable,
129 // with weight (by default w=1)
130 //
131 fData->Fill(var,weight);
132}
133
134//___________________________________________________________________
135TH1D *AliCFGridSparse::Project(Int_t ivar) const
136{
137 //
138 // Make a 1D projection along variable ivar
139 //
140
7411edfd 141 //exclude overflows in hidden dimesniosn by default (used in projections)
318f64b1 142 if(fExclOffEntriesInProj){
db6722a5 143 for(Int_t i=0;i<fNVar;i++){
7411edfd 144 if(i!=ivar)
145 fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
db6722a5 146 }
147 }
148 TH1D *hist=fData->Projection(ivar);
db6722a5 149 return hist;
db6722a5 150}
151//___________________________________________________________________
152TH2D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2) const
153{
154 //
155 // Make a 2D projection along variables ivar1 & ivar2
156 //
157
158 //exclude overflows by default (used in projections)
318f64b1 159 if(fExclOffEntriesInProj){
db6722a5 160 for(Int_t i=0;i<fNVar;i++){
7411edfd 161 if(i!=ivar1 && i!=ivar2)
162 fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
db6722a5 163 }
164 }
318f64b1 165 TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
db6722a5 166 return hist;
167
168}
169//___________________________________________________________________
170TH3D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
171{
172 //
173 // Make a 3D projection along variables ivar1 & ivar2 & ivar3
174 //
175 //exclude overflows by default (used in projections)
318f64b1 176 if(fExclOffEntriesInProj){
db6722a5 177 for(Int_t i=0;i<fNVar;i++){
7411edfd 178 if(i!=ivar1 && i!=ivar2 && i!=ivar3)
179 fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
db6722a5 180 }
181 }
182
183 TH3D *hist=fData->Projection(ivar1,ivar2,ivar3);
db6722a5 184 return hist;
185
186}
187
188//____________________________________________________________________
189Float_t AliCFGridSparse::GetOverFlows(Int_t ivar) const
190{
191 //
318f64b1 192 // Returns exclusive overflows in variable ivar
db6722a5 193 //
7411edfd 194 Int_t* bin = new Int_t[fNVar];
195 memset(bin, 0, sizeof(Int_t) * fNVar);
db6722a5 196 Float_t ovfl=0.;
197 for (Long64_t i = 0; i < fData->GetNbins(); ++i) {
198 Double_t v = fData->GetBinContent(i, bin);
199 Bool_t add=kTRUE;
200 for(Int_t j=0;j<fNVar;j++){
201 if(ivar==j)continue;
202 if((bin[j]==0) || (bin[j]==fNVarBins[j]+1))add=kFALSE;
203 }
204 if(bin[ivar]==fNVarBins[ivar]+1 && add) ovfl+=v;
205 }
206
207 delete[] bin;
208 return ovfl;
209}
210
211//____________________________________________________________________
212Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar) const
213{
214 //
318f64b1 215 // Returns exclusive overflows in variable ivar
db6722a5 216 //
7411edfd 217 Int_t* bin = new Int_t[fNVar];
218 memset(bin, 0, sizeof(Int_t) * fNVar);
db6722a5 219 Float_t unfl=0.;
220 for (Long64_t i = 0; i < fData->GetNbins(); ++i) {
221 Double_t v = fData->GetBinContent(i, bin);
222 Bool_t add=kTRUE;
223 for(Int_t j=0;j<fNVar;j++){
224 if(ivar==j)continue;
225 if((bin[j]==0) || (bin[j]==fNVarBins[j]+1))add=kFALSE;
226 }
227 if(bin[ivar]==0 && add) unfl+=v;
228 }
229
230 delete[] bin;
231 return unfl;
232}
233
234
235//____________________________________________________________________
236Float_t AliCFGridSparse::GetEntries() const
237{
238 //
239 // total entries (including overflows and underflows)
240 //
241
242 return fData->GetEntries();
243}
244
245//____________________________________________________________________
246Float_t AliCFGridSparse::GetElement(Int_t index) const
247{
248 //
249 // Returns content of grid element index according to the
250 // linear indexing in AliCFFrame
251 //
252 Int_t *bin = new Int_t[fNVar];
253 GetBinIndex(index, bin);
254 for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1; //consistency with AliCFGrid
255 Float_t val=GetElement(fIndex);
256 delete [] bin;
257 return val;
258
259}
260//____________________________________________________________________
261Float_t AliCFGridSparse::GetElement(Int_t *bin) const
262{
263 //
264 // Get the content in a bin corresponding to a set of bin indexes
265 //
266 return fData->GetBinContent(bin);
267
268}
269//____________________________________________________________________
270Float_t AliCFGridSparse::GetElement(Double_t *var) const
271{
272 //
273 // Get the content in a bin corresponding to a set of input variables
274 //
275
276 Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empty)
277 if(index<0){
278 return 0.;
279 }else{
280 return fData->GetBinContent(index);
281 }
282}
283
284//____________________________________________________________________
285Float_t AliCFGridSparse::GetElementError(Int_t index) const
286{
287 //
288 // Returns the error on the content of a bin according to a linear
289 // indexing in AliCFFrame
290 //
291
292 Int_t *bin = new Int_t[fNVar];
293 GetBinIndex(index, bin);
294 for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1; //consistency with AliCFGrid
295 Float_t val=GetElementError(fIndex);
296 delete [] bin;
297 return val;
298
299}
300//____________________________________________________________________
301Float_t AliCFGridSparse::GetElementError(Int_t *bin) const
302{
303 //
304 // Get the error in a bin corresponding to a set of bin indexes
305 //
306 return fData->GetBinError(bin);
307
308}
309//____________________________________________________________________
310Float_t AliCFGridSparse::GetElementError(Double_t *var) const
311{
312 //
313 // Get the error in a bin corresponding to a set of input variables
314 //
315
316 Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empy)
317 if(index<0){
318 return 0.;
319 }else{
320 return fData->GetBinError(index);
321 }
322}
323
324
325//____________________________________________________________________
326void AliCFGridSparse::SetElement(Int_t index, Float_t val)
327{
328 //
329 // Sets grid element iel to val (linear indexing) in AliCFFrame
330 //
331 Int_t *bin = new Int_t[fNVar];
332 GetBinIndex(index, bin);
333 for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1;
334 SetElement(fIndex,val);
335 delete [] bin;
336}
337//____________________________________________________________________
338void AliCFGridSparse::SetElement(Int_t *bin, Float_t val)
339{
340 //
341 // Sets grid element of bin indeces bin to val
342 //
343 fData->SetBinContent(bin,val);
344}
345//____________________________________________________________________
346void AliCFGridSparse::SetElement(Double_t *var, Float_t val)
347{
348 //
349 // Set the content in a bin to value val corresponding to a set of input variables
350 //
351 Long_t index=fData->GetBin(var); //THnSparse index: allocate the cell
352 Int_t *bin = new Int_t[fNVar];
353 fData->GetBinContent(index,bin); //trick to access the array of bins
354 fData->SetBinContent(bin,val);
355 delete [] bin;
356
357}
358
359//____________________________________________________________________
360void AliCFGridSparse::SetElementError(Int_t index, Float_t val)
361{
362 //
363 // Sets grid element iel error to val (linear indexing) in AliCFFrame
364 //
365 Int_t *bin = new Int_t[fNVar];
366 GetBinIndex(index, bin);
367 for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1;
368 SetElementError(fIndex,val);
369 delete [] bin;
370}
371//____________________________________________________________________
372void AliCFGridSparse::SetElementError(Int_t *bin, Float_t val)
373{
374 //
375 // Sets grid element error of bin indeces bin to val
376 //
377 fData->SetBinError(bin,val);
378}
379//____________________________________________________________________
380void AliCFGridSparse::SetElementError(Double_t *var, Float_t val)
381{
382 //
383 // Set the error in a bin to value val corresponding to a set of input variables
384 //
385 Long_t index=fData->GetBin(var); //THnSparse index
386 Int_t *bin = new Int_t[fNVar];
387 fData->GetBinContent(index,bin); //trick to access the array of bins
388 fData->SetBinError(bin,val);
389 delete [] bin;
390}
391
392//____________________________________________________________________
393void AliCFGridSparse::SumW2()
394{
395 //
396 //set calculation of the squared sum of the weighted entries
397 //
398 if(!fSumW2){
399 fData->CalculateErrors(kTRUE);
400 }
401
402 fSumW2=kTRUE;
403}
404
405//____________________________________________________________________
406void AliCFGridSparse::Add(AliCFVGrid* aGrid, Double_t c)
407{
408 //
409 //add aGrid to the current one
410 //
411
412
413 if(aGrid->GetNVar()!=fNVar){
414 AliInfo("Different number of variables, cannot add the grids");
415 return;
416 }
417 if(aGrid->GetNDim()!=fNDim){
418 AliInfo("Different number of dimensions, cannot add the grids!");
419 return;
420 }
421
422 if(!fSumW2 && aGrid->GetSumW2())SumW2();
423
424
425 fData->Add(((AliCFGridSparse*)aGrid)->GetGrid(),c);
426
427}
428
429//____________________________________________________________________
430void AliCFGridSparse::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
431{
432 //
433 //Add aGrid1 and aGrid2 and deposit the result into the current one
434 //
435
436 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
437 AliInfo("Different number of variables, cannot add the grids");
438 return;
439 }
440 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
441 AliInfo("Different number of dimensions, cannot add the grids!");
442 return;
443 }
444
445 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
446
447
448 fData->Reset();
449 fData->Add(((AliCFGridSparse*)aGrid1)->GetGrid(),c1);
450 fData->Add(((AliCFGridSparse*)aGrid2)->GetGrid(),c2);
451
452}
453
454//____________________________________________________________________
455void AliCFGridSparse::Multiply(AliCFVGrid* aGrid, Double_t c)
456{
457 //
458 // Multiply aGrid to the current one
459 //
460
461
462 if(aGrid->GetNVar()!=fNVar){
463 AliInfo("Different number of variables, cannot multiply the grids");
464 return;
465 }
466 if(aGrid->GetNDim()!=fNDim){
467 AliInfo("Different number of dimensions, cannot multiply the grids!");
468 return;
469 }
470
471 if(!fSumW2 && aGrid->GetSumW2())SumW2();
472
473 THnSparse *h= ((AliCFGridSparse*)aGrid)->GetGrid();
474
475 fData->Multiply(h);
476 fData->Scale(c);
477
478}
479
480//____________________________________________________________________
481void AliCFGridSparse::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
482{
483 //
484 //Multiply aGrid1 and aGrid2 and deposit the result into the current one
485 //
486
487 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
488 AliInfo("Different number of variables, cannot multiply the grids");
489 return;
490 }
491 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
492 AliInfo("Different number of dimensions, cannot multiply the grids!");
493 return;
494 }
495
496 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
497
498
499 fData->Reset();
500 THnSparse *h1= ((AliCFGridSparse*)aGrid1)->GetGrid();
501 THnSparse *h2= ((AliCFGridSparse*)aGrid2)->GetGrid();
502 h2->Multiply(h1);
503 h2->Scale(c1*c2);
504 fData->Add(h2);
505}
506
507
508
509//____________________________________________________________________
510void AliCFGridSparse::Divide(AliCFVGrid* aGrid, Double_t c)
511{
512 //
513 // Divide aGrid to the current one
514 //
515
516
517 if(aGrid->GetNVar()!=fNVar){
518 AliInfo("Different number of variables, cannot divide the grids");
519 return;
520 }
521 if(aGrid->GetNDim()!=fNDim){
522 AliInfo("Different number of dimensions, cannot divide the grids!");
523 return;
524 }
525
526 if(!fSumW2 && aGrid->GetSumW2())SumW2();
527
528 THnSparse *h= ((AliCFGridSparse*)aGrid)->GetGrid();
529
530 fData->Divide(h);
531 fData->Scale(c);
532
533}
534
535//____________________________________________________________________
536void AliCFGridSparse::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
537{
538 //
539 //Divide aGrid1 and aGrid2 and deposit the result into the current one
540 //bynomial errors are supported
541 //
542
543 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
544 AliInfo("Different number of variables, cannot divide the grids");
545 return;
546 }
547 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
548 AliInfo("Different number of dimensions, cannot divide the grids!");
549 return;
550 }
551
552 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
553
554
555 THnSparse *h1= ((AliCFGridSparse*)aGrid1)->GetGrid();
556 THnSparse *h2= ((AliCFGridSparse*)aGrid2)->GetGrid();
557 fData->Divide(h1,h2,c1,c2,option);
558}
559
560
561//____________________________________________________________________
7411edfd 562void AliCFGridSparse::Rebin(const Int_t* group)
563{
564 //
565 // rebin the grid according to Rebin() as in THnSparse
566 // Please notice that the original number of bins on
567 // a given axis has to be divisible by the rebin group.
568 //
569
570 for(Int_t i=0;i<fNVar;i++){
571 if(group[i]!=1)AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
572 }
573
574 THnSparse *rebinned =fData->Rebin(group);
575 fData->Reset();
576 fData = rebinned;
577
578 //redefine the needed stuff
579
580 Int_t ndimTot=1;
581 Int_t nbinTot=0;
582
583 //number of bins in each dimension, auxiliary variables
584
585 for(Int_t ivar=0;ivar<fNVar;ivar++){
586 Int_t nbins = fData->GetAxis(ivar)->GetNbins();
587 fNVarBins[ivar]=nbins;
588 ndimTot*=fNVarBins[ivar];
589 nbinTot+=(fNVarBins[ivar]+1);
590 Int_t offset=0;
591 for(Int_t i =0;i<ivar;i++)offset+=(fNVarBins[i]+1);
592 fOffset[ivar]=offset;
593 Int_t prod=1;
594 for(Int_t i=0;i<ivar;i++)prod*=fNVarBins[i];
595 fProduct[ivar]=prod;
596 }
597
598 fNDim=ndimTot;
599
600 //now the array of bin limits
601
602 delete fVarBinLimits;
603 fNVarBinLimits=nbinTot;
604 fVarBinLimits=new Double_t[fNVarBinLimits];
605
606 for(Int_t ivar=0;ivar<fNVar;ivar++){
607 Double_t low = fData->GetAxis(ivar)->GetXmin();
608 Double_t high = fData->GetAxis(ivar)->GetXmax();
609 const TArrayD *xbins = fData->GetAxis(ivar)->GetXbins();
610 if (xbins->fN == 0){
611 for(Int_t ibin=0;ibin<=fNVarBins[ivar];ibin++){
612 fVarBinLimits[ibin+fOffset[ivar]] = low + ibin*(high-low)/((Double_t) fNVarBins[ivar]);
613 }
614 }
615 else{
616
617 for(Int_t ibin=0;ibin<=fNVarBins[ivar];ibin++) {
618 fVarBinLimits[ibin+fOffset[ivar]] = xbins->At(ibin);
619 }
620 }
621 }
622
623}
624//____________________________________________________________________
db6722a5 625void AliCFGridSparse::Copy(TObject& c) const
626{
627 //
628 // copy function
629 //
630 AliCFGridSparse& target = (AliCFGridSparse &) c;
631
632 if(fData)target.fData = fData;
633}
634