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