Implemented AliCFContainer::ShowSlice and AliCFGridSparse::Slice
[u/mrichter/AliRoot.git] / CORRFW / AliCFGridSparse.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 // 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 //____________________________________________________________________
38 ClassImp(AliCFGridSparse)
39
40 //____________________________________________________________________
41 AliCFGridSparse::AliCFGridSparse() : 
42   AliCFVGrid(),
43   fData(0x0)
44 {
45   // default constructor
46 }
47 //____________________________________________________________________
48 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) : 
49   AliCFVGrid(name,title),
50   fData(0x0)
51 {
52   // default constructor
53 }
54 //____________________________________________________________________
55 AliCFGridSparse::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),
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 //____________________________________________________________________
78 AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) : 
79   AliCFVGrid(c),
80   fData(c.fData)
81 {
82   //
83   // copy constructor
84   //
85   ((AliCFGridSparse &)c).Copy(*this);
86 }
87
88 //____________________________________________________________________
89 AliCFGridSparse::~AliCFGridSparse()
90 {
91   //
92   // destructor
93   //
94   if(fData) delete fData;
95 }
96
97 //____________________________________________________________________
98 AliCFGridSparse &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 //____________________________________________________________________
109 void 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 //____________________________________________________________________
124 void 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 //___________________________________________________________________
135 TH1D *AliCFGridSparse::Project(Int_t ivar) const
136 {
137   //
138   // Make a 1D projection along variable ivar 
139   //
140
141   //exclude overflows in hidden dimesniosn by default (used in projections)
142   if(fExclOffEntriesInProj){
143     for(Int_t i=0;i<fNVar;i++){
144       if(i!=ivar)
145         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
146     }
147   }
148   TH1D *hist=fData->Projection(ivar);
149   return hist;
150 }
151 //___________________________________________________________________
152 TH2D *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)
159   if(fExclOffEntriesInProj){
160     for(Int_t i=0;i<fNVar;i++){
161       if(i!=ivar1 && i!=ivar2)
162         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
163     }
164   }
165   TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
166   return hist;
167
168 }
169 //___________________________________________________________________
170 TH3D *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)
176   if(fExclOffEntriesInProj){
177     for(Int_t i=0;i<fNVar;i++){
178       if(i!=ivar1 && i!=ivar2 && i!=ivar3)
179         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
180     }
181   }
182
183   TH3D *hist=fData->Projection(ivar1,ivar2,ivar3); 
184   return hist;
185
186 }
187
188 //____________________________________________________________________
189 Float_t AliCFGridSparse::GetOverFlows(Int_t ivar) const
190 {
191   //
192   // Returns exclusive overflows in variable ivar
193   //
194   Int_t* bin = new Int_t[fNVar];
195   memset(bin, 0, sizeof(Int_t) * fNVar);
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 //____________________________________________________________________
212 Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar) const
213 {
214   //
215   // Returns exclusive overflows in variable ivar
216   //
217   Int_t* bin = new Int_t[fNVar];
218   memset(bin, 0, sizeof(Int_t) * fNVar);
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 //____________________________________________________________________
236 Float_t AliCFGridSparse::GetEntries() const
237 {
238   //
239   // total entries (including overflows and underflows)
240   //
241
242   return fData->GetEntries();
243 }
244
245 //____________________________________________________________________
246 Float_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 //____________________________________________________________________
261 Float_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 //____________________________________________________________________
270 Float_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 //____________________________________________________________________
285 Float_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 //____________________________________________________________________
301 Float_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 //____________________________________________________________________
310 Float_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 //____________________________________________________________________
326 void 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 //____________________________________________________________________
338 void 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 //____________________________________________________________________
346 void 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 //____________________________________________________________________
360 void 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 //____________________________________________________________________
372 void 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 //____________________________________________________________________
380 void 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 //____________________________________________________________________
393 void 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 //____________________________________________________________________
406 void 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 //____________________________________________________________________
430 void 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 //____________________________________________________________________
455 void 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 //____________________________________________________________________
481 void 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 //____________________________________________________________________
510 void 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 //____________________________________________________________________
536 void 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 //____________________________________________________________________
562 void 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 //____________________________________________________________________
625 void 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
635 //____________________________________________________________________
636 TH1D* AliCFGridSparse::Slice(Int_t iVar, Double_t *varMin, Double_t *varMax) const
637 {
638   //
639   // return a slice (1D-projection) on variable iVar while axis ranges are defined with varMin,varMax
640   // arrays varMin and varMax contain the min and max values of each variable.
641   // therefore varMin and varMax must have their dimensions equal to fNVar
642   //
643   
644   THnSparse* clone = (THnSparse*)fData->Clone();
645   for (Int_t iAxis=0; iAxis<fNVar; iAxis++) {
646     clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
647   }
648   return clone->Projection(iVar);
649 }
650
651 //____________________________________________________________________
652 TH2D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Double_t *varMin, Double_t *varMax) const
653 {
654   //
655   // return a slice (2D-projection) on variables iVar1 and iVar2 while axis ranges are defined with varMin,varMax
656   // arrays varMin and varMax contain the min and max values of each variable.
657   // therefore varMin and varMax must have their dimensions equal to fNVar
658   //
659   
660   THnSparse* clone = (THnSparse*)fData->Clone();
661   for (Int_t iAxis=0; iAxis<fNVar; iAxis++) {
662     clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
663   }
664   return clone->Projection(iVar1,iVar2);
665 }
666
667 //____________________________________________________________________
668 TH3D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, Double_t *varMin, Double_t *varMax) const
669 {
670   //
671   // return a slice (3D-projection) on variables iVar1, iVar2 and iVar3 while axis ranges are defined with varMin,varMax
672   // arrays varMin and varMax contain the min and max values of each variable.
673   // therefore varMin and varMax must have their dimensions equal to fNVar
674   //
675
676   THnSparse* clone = (THnSparse*)fData->Clone();
677   for (Int_t iAxis=0; iAxis<fNVar; iAxis++) {
678     clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
679   }
680   return clone->Projection(iVar1,iVar2,iVar3);
681 }
682
683 //____________________________________________________________________
684 void AliCFGridSparse::SetRangeUser(Double_t *varMin, Double_t *varMax) {
685   //
686   // set range of every axis. varMin and varMax must be of dimension fNVar
687   //
688   for (Int_t iAxis=0; iAxis<fNVar ; iAxis++) { // set new range for every axis
689     fData->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
690   }
691   AliWarning("THnSparse has been modified");
692 }