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