1c4f31232032362146a26e7e4fd4fe5a2d299679
[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   AliCFFrame(),
43   fSumW2(kFALSE),
44   fData(0x0)
45 {
46   // default constructor
47 }
48 //____________________________________________________________________
49 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) : 
50   AliCFFrame(name,title),
51   fSumW2(kFALSE),
52   fData(0x0)
53 {
54   // default constructor
55 }
56 //____________________________________________________________________
57 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, Int_t nVarIn, const Int_t * nBinIn) :  
58   AliCFFrame(name,title),
59   fSumW2(kFALSE),
60   fData(0x0)
61 {
62   //
63   // main constructor
64   //
65
66   fData = new THnSparseF(name,title,nVarIn,nBinIn);
67 }
68
69 //____________________________________________________________________
70 AliCFGridSparse::~AliCFGridSparse()
71 {
72   //
73   // destructor
74   //
75   if (fData) delete fData;
76 }
77
78 //____________________________________________________________________
79 AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
80   AliCFFrame(c),
81   fSumW2(kFALSE),
82   fData(0x0)
83 {
84   //
85   // copy constructor
86   //
87   ((AliCFGridSparse &)c).Copy(*this);
88 }
89
90 //____________________________________________________________________
91 AliCFGridSparse& AliCFGridSparse::operator=(const AliCFGridSparse &c)
92 {
93   //
94   // assigment operator
95   //
96   if (this != &c) c.Copy(*this);
97   return *this;
98
99
100 //____________________________________________________________________
101 void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t min, Double_t max)
102 {
103   //
104   // set a uniform binning for variable ivar
105   //
106   Int_t nBins = GetNBins(ivar);
107   Double_t * array = new Double_t[nBins+1];
108   for (Int_t iEdge=0; iEdge<=nBins; iEdge++) array[iEdge] = min + iEdge * (max-min)/nBins ;
109   fData->SetBinEdges(ivar, array);
110   delete [] array ;
111
112
113 //____________________________________________________________________
114 void AliCFGridSparse::SetBinLimits(Int_t ivar, const Double_t *array)
115 {
116   //
117   // setting the arrays containing the bin limits 
118   //
119   fData->SetBinEdges(ivar, array);
120
121
122 //____________________________________________________________________
123 void AliCFGridSparse::Fill(const Double_t *var, Double_t weight)
124 {
125   //
126   // Fill the grid,
127   // given a set of values of the input variable, 
128   // with weight (by default w=1)
129   //
130   fData->Fill(var,weight);
131 }
132
133 //___________________________________________________________________
134 TH1D *AliCFGridSparse::Project(Int_t ivar) const
135 {
136   //
137   // Make a 1D projection along variable ivar 
138   //
139
140   TH1D *hist=fData->Projection(ivar);
141   hist->SetXTitle(GetVarTitle(ivar));
142   hist->SetName(Form("%s_proj-%s",GetName(),GetVarTitle(ivar)));
143   hist->SetTitle(Form("%s: projection on %s",GetTitle(),GetVarTitle(ivar)));
144   for (Int_t iBin=1; iBin<=GetNBins(ivar); iBin++) {
145     TString binLabel = GetAxis(ivar)->GetBinLabel(iBin) ;
146     if (binLabel.CompareTo("") != 0) hist->GetXaxis()->SetBinLabel(iBin,binLabel);
147   }
148   return hist;
149 }
150 //___________________________________________________________________
151 TH2D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2) const
152 {
153   //
154   // Make a 2D projection along variables ivar1 & ivar2 
155   //
156
157   TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
158   hist->SetXTitle(GetVarTitle(ivar1));
159   hist->SetYTitle(GetVarTitle(ivar2));
160   hist->SetName(Form("%s_proj-%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
161   hist->SetTitle(Form("%s: projection on %s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
162   for (Int_t iBin=1; iBin<=GetNBins(ivar1); iBin++) {
163     TString binLabel = GetAxis(ivar1)->GetBinLabel(iBin) ;
164     if (binLabel.CompareTo("") != 0) hist->GetXaxis()->SetBinLabel(iBin,binLabel);
165   }
166   for (Int_t iBin=1; iBin<=GetNBins(ivar2); iBin++) {
167     TString binLabel = GetAxis(ivar2)->GetBinLabel(iBin) ;
168     if (binLabel.CompareTo("") != 0) hist->GetYaxis()->SetBinLabel(iBin,binLabel);
169   }
170   return hist;
171 }
172 //___________________________________________________________________
173 TH3D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
174 {
175   //
176   // Make a 3D projection along variables ivar1 & ivar2 & ivar3 
177   //
178
179   TH3D *hist=fData->Projection(ivar1,ivar2,ivar3); 
180   hist->SetXTitle(GetVarTitle(ivar1));
181   hist->SetYTitle(GetVarTitle(ivar2));
182   hist->SetZTitle(GetVarTitle(ivar3));
183   hist->SetName(Form("%s_proj-%s,%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
184   hist->SetTitle(Form("%s: projection on %s-%s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
185   for (Int_t iBin=1; iBin<=GetNBins(ivar1); iBin++) {
186     TString binLabel = GetAxis(ivar1)->GetBinLabel(iBin) ;
187     if (binLabel.CompareTo("") != 0) hist->GetXaxis()->SetBinLabel(iBin,binLabel);
188   }
189   for (Int_t iBin=1; iBin<=GetNBins(ivar2); iBin++) {
190     TString binLabel = GetAxis(ivar2)->GetBinLabel(iBin) ;
191     if (binLabel.CompareTo("") != 0) hist->GetYaxis()->SetBinLabel(iBin,binLabel);
192   }
193   for (Int_t iBin=1; iBin<=GetNBins(ivar3); iBin++) {
194     TString binLabel = GetAxis(ivar3)->GetBinLabel(iBin) ;
195     if (binLabel.CompareTo("") != 0) hist->GetZaxis()->SetBinLabel(iBin,binLabel);
196   }
197   return hist;
198 }
199
200 //___________________________________________________________________
201 AliCFGridSparse* AliCFGridSparse::Project(Int_t nVars, const Int_t* vars, const Double_t* varMin, const Double_t* varMax, Bool_t useBins) const
202 {
203   //
204   // projects the grid on the nVars dimensions defined in vars.
205   // axis ranges can be defined in arrays varMin, varMax
206   // If useBins=true, varMin and varMax are taken as bin numbers
207   //
208
209   // binning for new grid
210   Int_t* bins = new Int_t[nVars];
211   for (Int_t iVar=0; iVar<nVars; iVar++) {
212     bins[iVar] = GetNBins(vars[iVar]);
213   }
214   
215   // create new grid sparse
216   AliCFGridSparse* out = new AliCFGridSparse(fName,fTitle,nVars,bins);
217
218   //set the range in the THnSparse to project
219   THnSparse* clone = ((THnSparse*)fData->Clone());
220   if (varMin && varMax) {
221     for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
222       if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
223       else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
224     }
225   }
226   else AliInfo("Keeping same axis ranges");
227
228   out->SetGrid(clone->Projection(nVars,vars));
229   delete [] bins;
230   return out;
231 }
232
233
234 //____________________________________________________________________
235 Float_t AliCFGridSparse::GetBinCenter(Int_t ivar, Int_t ibin) const
236 {
237   //
238   // Returns the center of specified bin for variable axis ivar
239   // 
240   
241   return (Float_t) fData->GetAxis(ivar)->GetBinCenter(ibin);
242 }
243
244 //____________________________________________________________________
245 Float_t AliCFGridSparse::GetBinSize(Int_t ivar, Int_t ibin) const
246 {
247   //
248   // Returns the size of specified bin for variable axis ivar
249   // 
250   
251   return (Float_t) fData->GetAxis(ivar)->GetBinUpEdge(ibin) - fData->GetAxis(ivar)->GetBinLowEdge(ibin);
252 }
253
254 //____________________________________________________________________
255 Float_t AliCFGridSparse::GetEntries() const
256 {
257   //
258   // total entries (including overflows and underflows)
259   //
260
261   return fData->GetEntries();
262 }
263
264 //____________________________________________________________________
265 Float_t AliCFGridSparse::GetElement(Long_t index) const
266 {
267   //
268   // Returns content of grid element index 
269   //
270   
271   return fData->GetBinContent(index);
272 }
273 //____________________________________________________________________
274 Float_t AliCFGridSparse::GetElement(const Int_t *bin) const
275 {
276   //
277   // Get the content in a bin corresponding to a set of bin indexes
278   //
279   return fData->GetBinContent(bin);
280
281 }  
282 //____________________________________________________________________
283 Float_t AliCFGridSparse::GetElement(const Double_t *var) const
284 {
285   //
286   // Get the content in a bin corresponding to a set of input variables
287   //
288
289   Long_t index = fData->GetBin(var,kFALSE);
290   if (index<0) return 0.;
291   return fData->GetBinContent(index);
292
293
294 //____________________________________________________________________
295 Float_t AliCFGridSparse::GetElementError(Long_t index) const
296 {
297   //
298   // Returns the error on the content 
299   //
300
301   return fData->GetBinContent(index);
302 }
303 //____________________________________________________________________
304 Float_t AliCFGridSparse::GetElementError(const Int_t *bin) const
305 {
306  //
307   // Get the error in a bin corresponding to a set of bin indexes
308   //
309   return fData->GetBinError(bin);
310
311 }  
312 //____________________________________________________________________
313 Float_t AliCFGridSparse::GetElementError(const Double_t *var) const
314 {
315   //
316   // Get the error in a bin corresponding to a set of input variables
317   //
318
319   Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empy)
320   if (index<0) return 0.;
321   return fData->GetBinError(index);
322
323
324 //____________________________________________________________________
325 void AliCFGridSparse::SetElement(Long_t index, Float_t val)
326 {
327   //
328   // Sets grid element value
329   //
330   Int_t* bin = new Int_t[GetNVar()];
331   fData->GetBinContent(index,bin); //affects the bin coordinates
332   SetElement(bin,val);
333   delete [] bin ;
334 }
335
336 //____________________________________________________________________
337 void AliCFGridSparse::SetElement(const Int_t *bin, Float_t val)
338 {
339   //
340   // Sets grid element of bin indeces bin to val
341   //
342   fData->SetBinContent(bin,val);
343 }
344 //____________________________________________________________________
345 void AliCFGridSparse::SetElement(const Double_t *var, Float_t val) 
346 {
347   //
348   // Set the content in a bin to value val corresponding to a set of input variables
349   //
350   Long_t index=fData->GetBin(var,kTRUE); //THnSparse index: allocate the cell
351   Int_t *bin = new Int_t[GetNVar()];
352   fData->GetBinContent(index,bin); //trick to access the array of bins
353   SetElement(bin,val);
354   delete [] bin;
355 }
356
357 //____________________________________________________________________
358 void AliCFGridSparse::SetElementError(Long_t index, Float_t val)
359 {
360   //
361   // Sets grid element iel error to val (linear indexing) in AliCFFrame
362   //
363   Int_t *bin = new Int_t[GetNVar()];
364   fData->GetBinContent(index,bin);
365   SetElementError(bin,val);
366   delete [] bin;
367 }
368
369 //____________________________________________________________________
370 void AliCFGridSparse::SetElementError(const Int_t *bin, Float_t val)
371 {
372   //
373   // Sets grid element error of bin indeces bin to val
374   //
375   fData->SetBinError(bin,val);
376 }
377 //____________________________________________________________________
378 void AliCFGridSparse::SetElementError(const Double_t *var, Float_t val) 
379 {
380   //
381   // Set the error in a bin to value val corresponding to a set of input variables
382   //
383   Long_t index=fData->GetBin(var); //THnSparse index
384   Int_t *bin = new Int_t[GetNVar()];
385   fData->GetBinContent(index,bin); //trick to access the array of bins
386   SetElementError(bin,val);
387   delete [] bin;
388 }
389
390 //____________________________________________________________________
391 void AliCFGridSparse::SumW2()
392 {
393   //
394   //set calculation of the squared sum of the weighted entries
395   //
396   if(!fSumW2){
397     fData->CalculateErrors(kTRUE); 
398   }
399   fSumW2=kTRUE;
400 }
401
402 //____________________________________________________________________
403 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid, Double_t c)
404 {
405   //
406   //add aGrid to the current one
407   //
408
409   if (aGrid->GetNVar() != GetNVar()){
410     AliError("Different number of variables, cannot add the grids");
411     return;
412   } 
413   
414   if (!fSumW2  && aGrid->GetSumW2()) SumW2();
415   fData->Add(aGrid->GetGrid(),c);
416 }
417
418 //____________________________________________________________________
419 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
420 {
421   //
422   //Add aGrid1 and aGrid2 and deposit the result into the current one
423   //
424
425   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
426     AliInfo("Different number of variables, cannot add the grids");
427     return;
428   } 
429   
430   if (!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
431
432   fData->Reset();
433   fData->Add(aGrid1->GetGrid(),c1);
434   fData->Add(aGrid2->GetGrid(),c2);
435 }
436
437 //____________________________________________________________________
438 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid, Double_t c)
439 {
440   //
441   // Multiply aGrid to the current one
442   //
443
444   if (aGrid->GetNVar() != GetNVar()) {
445     AliError("Different number of variables, cannot multiply the grids");
446     return;
447   } 
448   
449   if(!fSumW2  && aGrid->GetSumW2()) SumW2();
450   THnSparse *h = aGrid->GetGrid();
451   fData->Multiply(h);
452   fData->Scale(c);
453 }
454
455 //____________________________________________________________________
456 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
457 {
458   //
459   //Multiply aGrid1 and aGrid2 and deposit the result into the current one
460   //
461
462   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
463     AliError("Different number of variables, cannot multiply the grids");
464     return;
465   }
466   
467   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
468
469   fData->Reset();
470   THnSparse *h1 = aGrid1->GetGrid();
471   THnSparse *h2 = aGrid2->GetGrid();
472   h2->Multiply(h1);
473   h2->Scale(c1*c2);
474   fData->Add(h2);
475 }
476
477 //____________________________________________________________________
478 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid, Double_t c)
479 {
480   //
481   // Divide aGrid to the current one
482   //
483
484   if (aGrid->GetNVar() != GetNVar()) {
485     AliError("Different number of variables, cannot divide the grids");
486     return;
487   } 
488   
489   if (!fSumW2  && aGrid->GetSumW2()) SumW2();
490
491   THnSparse *h1 = aGrid->GetGrid();
492   THnSparse *h2 = (THnSparse*)fData->Clone();
493   fData->Divide(h2,h1);
494   fData->Scale(c);
495 }
496
497 //____________________________________________________________________
498 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2, Option_t *option)
499 {
500   //
501   //Divide aGrid1 and aGrid2 and deposit the result into the current one
502   //binomial errors are supported
503   //
504
505   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
506     AliError("Different number of variables, cannot divide the grids");
507     return;
508   } 
509   
510   if (!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
511
512   THnSparse *h1= aGrid1->GetGrid();
513   THnSparse *h2= aGrid2->GetGrid();
514   fData->Divide(h1,h2,c1,c2,option);
515 }
516
517
518 //____________________________________________________________________
519 void AliCFGridSparse::Rebin(const Int_t* group)
520 {
521   //
522   // rebin the grid according to Rebin() as in THnSparse
523   // Please notice that the original number of bins on
524   // a given axis has to be divisible by the rebin group.
525   //
526
527   for(Int_t i=0;i<GetNVar();i++){
528     if (group[i]!=1) AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
529   }
530
531   THnSparse *rebinned =fData->Rebin(group);
532   fData->Reset();
533   fData = rebinned;
534 }
535 //____________________________________________________________________
536 void AliCFGridSparse::Scale(Long_t index, const Double_t *fact)
537 {
538   //
539   //scale content of a certain cell by (positive) fact (with error)
540   //
541
542   if (GetElement(index)==0 || fact[0]==0) return;
543
544   Double_t in[2], out[2];
545   in[0]=GetElement(index);
546   in[1]=GetElementError(index);
547   GetScaledValues(fact,in,out);
548   SetElement(index,out[0]);
549   if (fSumW2) SetElementError(index,out[1]);
550 }
551 //____________________________________________________________________
552 void AliCFGridSparse::Scale(const Int_t *bin, const Double_t *fact)
553 {
554   //
555   //scale content of a certain cell by (positive) fact (with error)
556   //
557   if(GetElement(bin)==0 || fact[0]==0)return;
558
559   Double_t in[2], out[2];
560   in[0]=GetElement(bin);
561   in[1]=GetElementError(bin);
562   GetScaledValues(fact,in,out);
563   SetElement(bin,out[0]);
564   if(fSumW2)SetElementError(bin,out[1]);
565   
566 }
567 //____________________________________________________________________
568 void AliCFGridSparse::Scale(const Double_t *var, const Double_t *fact) 
569 {
570   //
571   //scale content of a certain cell by (positive) fact (with error)
572   //
573   if(GetElement(var)==0 || fact[0]==0)return;
574
575   Double_t in[2], out[2];
576   in[0]=GetElement(var);
577   in[1]=GetElementError(var);
578   GetScaledValues(fact,in,out);
579   SetElement(var,out[0]);
580   if(fSumW2)SetElementError(var,out[1]);
581   
582 }
583 //____________________________________________________________________
584 void AliCFGridSparse::Scale(const Double_t* fact)
585 {
586   //
587   //scale contents of the whole grid by fact
588   //
589
590   for (Long_t iel=0; iel<GetNFilledBins(); iel++) {
591     Scale(iel,fact);
592   }
593 }
594 //____________________________________________________________________
595 Long_t AliCFGridSparse::GetEmptyBins() const {
596   //
597   // Get empty bins 
598   //
599
600   return (GetNBinsTotal() - GetNFilledBins()) ;
601
602
603 //_____________________________________________________________________
604 // Int_t AliCFGridSparse::GetEmptyBins( Double_t *varMin, Double_t* varMax ) const 
605 // {
606 //   //
607 //   // Get empty bins in a range specified by varMin and varMax
608 //   //
609
610 //   Int_t *indexMin = new Int_t[GetNVar()];
611 //   Int_t *indexMax = new Int_t[GetNVar()];
612
613 //   //Find out the min and max bins
614
615 //   for (Int_t i=0;i<GetNVar();i++) {
616 //     Double_t xmin=varMin[i]; // the min values  
617 //     Double_t xmax=varMax[i]; // the min values  
618 //     Int_t nbins = fData->GetAxis(i)->GetNbins()+1;
619 //     Double_t *bins=new Double_t[nbins];
620 //     for(Int_t ibin =0;ibin<nbins;ibin++){
621 //      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
622 //     }
623 //     indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
624 //     indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
625 //     delete [] bins;
626 //   }
627
628 //   Int_t val=0;
629 //   for(Int_t i=0;i<fNDim;i++){
630 //     for (Int_t j=0;j<GetNVar();j++)fIndex[j]=GetBinIndex(j,i);
631 //     Bool_t isIn=kTRUE;
632 //     for (Int_t j=0;j<GetNVar();j++){
633 //       if(!(fIndex[j]>=indexMin[j] && fIndex[j]<=indexMax[j]))isIn=kFALSE;   
634 //     }
635 //     if(isIn && GetElement(i)<=0)val++;     
636 //   }
637 //   AliInfo(Form(" the empty bins = %i ",val)); 
638
639 //   delete [] indexMin;
640 //   delete [] indexMax;
641 //   return val;
642 // } 
643
644 //____________________________________________________________________
645 Int_t AliCFGridSparse::CheckStats(Double_t thr) const
646 {
647   //
648   // Count the cells below a certain threshold
649   //
650   Int_t ncellsLow=0;
651   for (Int_t i=0; i<GetNBinsTotal(); i++) {
652     if (GetElement(i)<thr) ncellsLow++;
653   }
654   return ncellsLow;
655 }
656
657 //_____________________________________________________________________
658 Double_t AliCFGridSparse::GetIntegral() const 
659 {
660   //
661   // Get full Integral
662   //
663   return fData->ComputeIntegral();  
664
665
666 //_____________________________________________________________________
667 // Double_t AliCFGridSparse::GetIntegral(Int_t *binMin, Int_t* binMax ) const 
668 // {
669 //   //
670 //   // Get Integral in a range of bin indeces (extremes included)
671 //   //
672
673 //   Double_t val=0;
674
675 //   for(Int_t i=0;i<GetNVar();i++){
676 //     if(binMin[i]<1)binMin[i]=1;
677 //     if(binMax[i]>fNVarBins[i])binMax[i]=fNVarBins[i];
678 //     if((binMin[i]>binMax[i])){
679 //       AliInfo(Form(" Bin indices in variable %i in reverse order, please check!", i));
680 //       return val;
681 //     }
682 //   }
683 //   val=GetSum(0,binMin,binMax);
684
685 //   return val;
686 // } 
687
688 //_____________________________________________________________________
689 // Double_t AliCFGridSparse::GetIntegral(Double_t *varMin, Double_t* varMax ) const 
690 // {
691 //   //
692 //   // Get Integral in a range (extremes included)
693 //   //
694
695 //   Int_t *indexMin=new Int_t[GetNVar()];
696 //   Int_t *indexMax=new Int_t[GetNVar()];
697
698 //   //Find out the min and max bins
699
700 //   for(Int_t i=0;i<GetNVar();i++){
701 //     Double_t xmin=varMin[i]; // the min values  
702 //     Double_t xmax=varMax[i]; // the min values  
703 //     Int_t nbins=fNVarBins[i]+1;
704 //     Double_t *bins=new Double_t[nbins];
705 //     for(Int_t ibin =0;ibin<nbins;ibin++){
706 //      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
707 //     }
708 //     indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
709 //     indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
710 //     delete [] bins;
711 //   }
712
713 //   //move to the TH/THnSparse convention in N-dim bin numbering 
714 //   for(Int_t i=0;i<GetNVar();i++){
715 //     indexMin[i]+=1;
716 //     indexMax[i]+=1;
717 //   }
718
719 //   Double_t val=GetIntegral(indexMin,indexMax);
720
721 //   delete [] indexMin;
722 //   delete [] indexMax;
723
724 //   return val;
725 // } 
726
727 // //_____________________________________________________________________
728 // Double_t AliCFGridSparse::GetSum(Int_t ivar, Int_t *binMin, Int_t* binMax) const 
729 // {
730 //   //
731 //   // recursively add over nested loops.... 
732 //   //
733
734 //   static Double_t val;
735 //   if (ivar==0) val=0.;
736 //   for(Int_t ibin=binMin[ivar]-1 ; ibin<=binMax[ivar]-1 ; ibin++){
737 //     //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
738 //     fIndex[ivar]=ibin;
739 //     if(ivar<GetNVar()-1) {
740 //       val=GetSum(ivar+1,binMin,binMax);
741 //     }
742 //     else {
743 //       Int_t iel=GetBinIndex(fIndex);
744 //       val+=GetElement(iel);
745 //     }
746 //   }
747   
748 //   return val;
749 // }
750
751 //____________________________________________________________________
752 Long64_t AliCFGridSparse::Merge(TCollection* list)
753 {
754   //
755   // Merge a list of AliCFGridSparse with this (needed for PROOF). 
756   // Returns the number of merged objects (including this).
757   //
758
759   if (!list)
760     return 0;
761   
762   if (list->IsEmpty())
763     return 1;
764
765   TIterator* iter = list->MakeIterator();
766   TObject* obj;
767   
768   Int_t count = 0;
769   while ((obj = iter->Next())) {
770     AliCFGridSparse* entry = dynamic_cast<AliCFGridSparse*> (obj);
771     if (entry == 0) 
772       continue;
773     this->Add(entry);
774     count++;
775   }
776
777   return count+1;
778 }
779
780 //____________________________________________________________________
781 void AliCFGridSparse::GetScaledValues(const Double_t *fact, const Double_t *in, Double_t *out) const{
782   //
783   // scale input *in and its error by (positive) fact (with error)
784   // and erite it to *out
785   //
786   out[0]=in[0]*fact[0];
787   out[1]=TMath::Sqrt(in[1]*in[1]/in[0]/in[0]
788                      +fact[1]*fact[1]/fact[0]/fact[0])*out[0];
789     
790 }
791
792 //____________________________________________________________________
793 void AliCFGridSparse::Copy(TObject& c) const
794 {
795   //
796   // copy function
797   //
798   AliCFFrame::Copy(c);
799   AliCFGridSparse& target = (AliCFGridSparse &) c;
800   target.fSumW2 = fSumW2 ;
801   if (fData) {
802     target.fData = (THnSparse*)fData->Clone();
803   }
804 }
805
806 //____________________________________________________________________
807 TH1D* AliCFGridSparse::Slice(Int_t iVar, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
808 {
809   //
810   // return a slice (1D-projection) on variable iVar while axis ranges are defined with varMin,varMax
811   // arrays varMin and varMax contain the min and max values of each variable.
812   // therefore varMin and varMax must have their dimensions equal to GetNVar()
813   // If useBins=true, varMin and varMax are taken as bin numbers
814   
815   THnSparse* clone = (THnSparse*)fData->Clone();
816   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
817     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
818     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
819   }
820
821   TH1D* projection = 0x0 ;
822
823   TObject* objInMem = gROOT->FindObject(Form("%s_proj_%d",clone->GetName(),iVar)) ;
824   if (objInMem) {
825     TString name(objInMem->ClassName());
826     if (name.CompareTo("TH1D")==0) projection = (TH1D*) objInMem ;
827     else projection = clone->Projection(iVar);
828   }
829   else projection = clone->Projection(iVar); 
830   delete clone;
831   for (Int_t iBin=1; iBin<=GetNBins(iVar); iBin++) {
832     TString binLabel = GetAxis(iVar)->GetBinLabel(iBin) ;
833     if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
834   }
835   return projection ;
836 }
837
838 //____________________________________________________________________
839 TH2D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
840 {
841   //
842   // return a slice (2D-projection) on variables iVar1 and iVar2 while axis ranges are defined with varMin,varMax
843   // arrays varMin and varMax contain the min and max values of each variable.
844   // therefore varMin and varMax must have their dimensions equal to GetNVar()
845   // If useBins=true, varMin and varMax are taken as bin numbers
846   
847   THnSparse* clone = (THnSparse*)fData->Clone();
848   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
849     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
850     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
851   }
852
853   TH2D* projection = 0x0 ;
854
855   TObject* objInMem = gROOT->FindObject(Form("%s_proj_%d_%d",clone->GetName(),iVar2,iVar1)) ;
856   if (objInMem) {
857     TString name(objInMem->ClassName());
858     if (name.CompareTo("TH2D")==0) projection = (TH2D*) objInMem ;
859     else projection = clone->Projection(iVar1,iVar2);
860   }
861   else projection = clone->Projection(iVar1,iVar2); 
862   delete clone;
863   for (Int_t iBin=1; iBin<=GetNBins(iVar1); iBin++) {
864     TString binLabel = GetAxis(iVar1)->GetBinLabel(iBin) ;
865     if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
866   }
867   for (Int_t iBin=1; iBin<=GetNBins(iVar2); iBin++) {
868     TString binLabel = GetAxis(iVar2)->GetBinLabel(iBin) ;
869     if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
870   }
871   return projection ;
872 }
873
874 //____________________________________________________________________
875 TH3D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
876 {
877   //
878   // return a slice (3D-projection) on variables iVar1, iVar2 and iVar3 while axis ranges are defined with varMin,varMax
879   // arrays varMin and varMax contain the min and max values of each variable.
880   // therefore varMin and varMax must have their dimensions equal to GetNVar()
881   // If useBins=true, varMin and varMax are taken as bin numbers
882
883   THnSparse* clone = (THnSparse*)fData->Clone();
884   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
885     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
886     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
887   }
888
889   TH3D* projection = 0x0 ;
890
891   TObject* objInMem = gROOT->FindObject(Form("%s_proj_%d_%d_%d",clone->GetName(),iVar1,iVar2,iVar3)) ;
892   if (objInMem) {
893     TString name(objInMem->ClassName());
894     if (name.CompareTo("TH3D")==0) projection = (TH3D*) objInMem ;
895     else projection = clone->Projection(iVar1,iVar2,iVar3);
896   }
897   else projection = clone->Projection(iVar1,iVar2,iVar3); 
898   delete clone;
899   for (Int_t iBin=1; iBin<=GetNBins(iVar1); iBin++) {
900     TString binLabel = GetAxis(iVar1)->GetBinLabel(iBin) ;
901     if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
902   }
903   for (Int_t iBin=1; iBin<=GetNBins(iVar2); iBin++) {
904     TString binLabel = GetAxis(iVar2)->GetBinLabel(iBin) ;
905     if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
906   }
907   for (Int_t iBin=1; iBin<=GetNBins(iVar3); iBin++) {
908     TString binLabel = GetAxis(iVar3)->GetBinLabel(iBin) ;
909     if (binLabel.CompareTo("") != 0) projection->GetZaxis()->SetBinLabel(iBin,binLabel);
910   }
911   return projection ;
912 }
913
914 //____________________________________________________________________
915 void AliCFGridSparse::SetRangeUser(Int_t iVar, Double_t varMin, Double_t varMax) {
916   //
917   // set range of axis iVar. 
918   //
919   fData->GetAxis(iVar)->SetRangeUser(varMin,varMax);
920   AliWarning(Form("THnSparse axis %d range has been modified",iVar));
921 }
922
923 //____________________________________________________________________
924 void AliCFGridSparse::SetRangeUser(const Double_t *varMin, const Double_t *varMax) {
925   //
926   // set range of every axis. varMin and varMax must be of dimension GetNVar()
927   //
928   for (Int_t iAxis=0; iAxis<GetNVar() ; iAxis++) { // set new range for every axis
929     SetRangeUser(iAxis,varMin[iAxis],varMax[iAxis]);
930   }
931   AliWarning("THnSparse axes ranges have been modified");
932 }
933
934 //____________________________________________________________________
935 void AliCFGridSparse::UseAxisRange(Bool_t b) const {
936   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) fData->GetAxis(iAxis)->SetBit(TAxis::kAxisRange,b);
937 }
938
939 //____________________________________________________________________
940 Float_t AliCFGridSparse::GetOverFlows(Int_t ivar, Bool_t exclusive) const
941 {
942   //
943   // Returns overflows in variable ivar
944   // Set 'exclusive' to true for an exclusive check on variable ivar
945   //
946   Int_t* bin = new Int_t[GetNVar()];
947   memset(bin, 0, sizeof(Int_t) * GetNVar());
948   Float_t ovfl=0.;
949   for (Long64_t i = 0; i < fData->GetNbins(); i++) {
950     Double_t v = fData->GetBinContent(i, bin);
951     Bool_t add=kTRUE;
952     if (exclusive) {
953       for(Int_t j=0;j<GetNVar();j++){
954         if(ivar==j)continue;
955         if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
956       }
957     }
958     if(bin[ivar]==GetNBins(ivar)+1 && add) ovfl+=v;
959   }
960
961   delete[] bin;
962   return ovfl;
963 }
964
965 //____________________________________________________________________
966 Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar, Bool_t exclusive) const
967 {
968   //
969   // Returns exclusive overflows in variable ivar
970   // Set 'exclusive' to true for an exclusive check on variable ivar
971   //
972   Int_t* bin = new Int_t[GetNVar()];
973   memset(bin, 0, sizeof(Int_t) * GetNVar());
974   Float_t unfl=0.;
975   for (Long64_t i = 0; i < fData->GetNbins(); i++) {
976     Double_t v = fData->GetBinContent(i, bin);
977     Bool_t add=kTRUE;
978     if (exclusive) {
979       for(Int_t j=0;j<GetNVar();j++){
980         if(ivar==j)continue;
981         if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
982       }
983     }
984     if(bin[ivar]==0 && add) unfl+=v;
985   }
986
987   delete[] bin;
988   return unfl;
989 }
990