1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 //--------------------------------------------------------------------//
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 //--------------------------------------------------------------------//
27 #include "AliCFGridSparse.h"
28 #include "THnSparse.h"
36 #include "AliCFUnfolding.h"
38 //____________________________________________________________________
39 ClassImp(AliCFGridSparse)
41 //____________________________________________________________________
42 AliCFGridSparse::AliCFGridSparse() :
47 // default constructor
49 //____________________________________________________________________
50 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
51 AliCFFrame(name,title),
55 // default constructor
57 //____________________________________________________________________
58 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, Int_t nVarIn, const Int_t * nBinIn) :
59 AliCFFrame(name,title),
67 fData = new THnSparseF(name,title,nVarIn,nBinIn);
70 //____________________________________________________________________
71 AliCFGridSparse::~AliCFGridSparse()
76 if (fData) delete fData;
79 //____________________________________________________________________
80 AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
88 ((AliCFGridSparse &)c).Copy(*this);
91 //____________________________________________________________________
92 AliCFGridSparse& AliCFGridSparse::operator=(const AliCFGridSparse &c)
97 if (this != &c) c.Copy(*this);
101 //____________________________________________________________________
102 void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t min, Double_t max)
105 // set a uniform binning for variable ivar
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);
114 //____________________________________________________________________
115 void AliCFGridSparse::SetBinLimits(Int_t ivar, const Double_t *array)
118 // setting the arrays containing the bin limits
120 fData->SetBinEdges(ivar, array);
123 //____________________________________________________________________
124 void AliCFGridSparse::Fill(const Double_t *var, Double_t weight)
128 // given a set of values of the input variable,
129 // with weight (by default w=1)
131 fData->Fill(var,weight);
134 //___________________________________________________________________
135 AliCFGridSparse* AliCFGridSparse::MakeSlice(Int_t nVars, const Int_t* vars, const Double_t* varMin, const Double_t* varMax, Bool_t useBins) const
138 // projects the grid on the nVars dimensions defined in vars.
139 // axis ranges can be defined in arrays varMin, varMax
140 // If useBins=true, varMin and varMax are taken as bin numbers
143 // binning for new grid
144 Int_t* bins = new Int_t[nVars];
145 for (Int_t iVar=0; iVar<nVars; iVar++) {
146 bins[iVar] = GetNBins(vars[iVar]);
149 // create new grid sparse
150 AliCFGridSparse* out = new AliCFGridSparse(fName,fTitle,nVars,bins);
152 //set the range in the THnSparse to project
153 THnSparse* clone = ((THnSparse*)fData->Clone());
154 if (varMin && varMax) {
155 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
156 SetAxisRange(clone->GetAxis(iAxis),varMin[iAxis],varMax[iAxis],useBins);
159 else AliInfo("Keeping same axis ranges");
161 out->SetGrid(clone->Projection(nVars,vars));
168 //____________________________________________________________________
169 Float_t AliCFGridSparse::GetBinCenter(Int_t ivar, Int_t ibin) const
172 // Returns the center of specified bin for variable axis ivar
175 return (Float_t) fData->GetAxis(ivar)->GetBinCenter(ibin);
178 //____________________________________________________________________
179 Float_t AliCFGridSparse::GetBinSize(Int_t ivar, Int_t ibin) const
182 // Returns the size of specified bin for variable axis ivar
185 return (Float_t) fData->GetAxis(ivar)->GetBinUpEdge(ibin) - fData->GetAxis(ivar)->GetBinLowEdge(ibin);
188 //____________________________________________________________________
189 Float_t AliCFGridSparse::GetEntries() const
192 // total entries (including overflows and underflows)
195 return fData->GetEntries();
198 //____________________________________________________________________
199 Float_t AliCFGridSparse::GetElement(Long_t index) const
202 // Returns content of grid element index
205 return fData->GetBinContent(index);
207 //____________________________________________________________________
208 Float_t AliCFGridSparse::GetElement(const Int_t *bin) const
211 // Get the content in a bin corresponding to a set of bin indexes
213 return fData->GetBinContent(bin);
216 //____________________________________________________________________
217 Float_t AliCFGridSparse::GetElement(const Double_t *var) const
220 // Get the content in a bin corresponding to a set of input variables
223 Long_t index = fData->GetBin(var,kFALSE);
224 if (index<0) return 0.;
225 return fData->GetBinContent(index);
228 //____________________________________________________________________
229 Float_t AliCFGridSparse::GetElementError(Long_t index) const
232 // Returns the error on the content
235 return fData->GetBinContent(index);
237 //____________________________________________________________________
238 Float_t AliCFGridSparse::GetElementError(const Int_t *bin) const
241 // Get the error in a bin corresponding to a set of bin indexes
243 return fData->GetBinError(bin);
246 //____________________________________________________________________
247 Float_t AliCFGridSparse::GetElementError(const Double_t *var) const
250 // Get the error in a bin corresponding to a set of input variables
253 Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empy)
254 if (index<0) return 0.;
255 return fData->GetBinError(index);
258 //____________________________________________________________________
259 void AliCFGridSparse::SetElement(Long_t index, Float_t val)
262 // Sets grid element value
264 Int_t* bin = new Int_t[GetNVar()];
265 fData->GetBinContent(index,bin); //affects the bin coordinates
270 //____________________________________________________________________
271 void AliCFGridSparse::SetElement(const Int_t *bin, Float_t val)
274 // Sets grid element of bin indeces bin to val
276 fData->SetBinContent(bin,val);
278 //____________________________________________________________________
279 void AliCFGridSparse::SetElement(const Double_t *var, Float_t val)
282 // Set the content in a bin to value val corresponding to a set of input variables
284 Long_t index=fData->GetBin(var,kTRUE); //THnSparse index: allocate the cell
285 Int_t *bin = new Int_t[GetNVar()];
286 fData->GetBinContent(index,bin); //trick to access the array of bins
291 //____________________________________________________________________
292 void AliCFGridSparse::SetElementError(Long_t index, Float_t val)
295 // Sets grid element iel error to val (linear indexing) in AliCFFrame
297 Int_t *bin = new Int_t[GetNVar()];
298 fData->GetBinContent(index,bin);
299 SetElementError(bin,val);
303 //____________________________________________________________________
304 void AliCFGridSparse::SetElementError(const Int_t *bin, Float_t val)
307 // Sets grid element error of bin indeces bin to val
309 fData->SetBinError(bin,val);
311 //____________________________________________________________________
312 void AliCFGridSparse::SetElementError(const Double_t *var, Float_t val)
315 // Set the error in a bin to value val corresponding to a set of input variables
317 Long_t index=fData->GetBin(var); //THnSparse index
318 Int_t *bin = new Int_t[GetNVar()];
319 fData->GetBinContent(index,bin); //trick to access the array of bins
320 SetElementError(bin,val);
324 //____________________________________________________________________
325 void AliCFGridSparse::SumW2()
328 //set calculation of the squared sum of the weighted entries
331 fData->CalculateErrors(kTRUE);
336 //____________________________________________________________________
337 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid, Double_t c)
340 //add aGrid to the current one
343 if (aGrid->GetNVar() != GetNVar()){
344 AliError("Different number of variables, cannot add the grids");
348 if (!fSumW2 && aGrid->GetSumW2()) SumW2();
349 fData->Add(aGrid->GetGrid(),c);
352 //____________________________________________________________________
353 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
356 //Add aGrid1 and aGrid2 and deposit the result into the current one
359 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
360 AliInfo("Different number of variables, cannot add the grids");
364 if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
367 fData->Add(aGrid1->GetGrid(),c1);
368 fData->Add(aGrid2->GetGrid(),c2);
371 //____________________________________________________________________
372 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid, Double_t c)
375 // Multiply aGrid to the current one
378 if (aGrid->GetNVar() != GetNVar()) {
379 AliError("Different number of variables, cannot multiply the grids");
383 if(!fSumW2 && aGrid->GetSumW2()) SumW2();
384 THnSparse *h = aGrid->GetGrid();
389 //____________________________________________________________________
390 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
393 //Multiply aGrid1 and aGrid2 and deposit the result into the current one
396 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
397 AliError("Different number of variables, cannot multiply the grids");
401 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
404 THnSparse *h1 = aGrid1->GetGrid();
405 THnSparse *h2 = aGrid2->GetGrid();
411 //____________________________________________________________________
412 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid, Double_t c)
415 // Divide aGrid to the current one
418 if (aGrid->GetNVar() != GetNVar()) {
419 AliError("Different number of variables, cannot divide the grids");
423 if (!fSumW2 && aGrid->GetSumW2()) SumW2();
425 THnSparse *h1 = aGrid->GetGrid();
426 THnSparse *h2 = (THnSparse*)fData->Clone();
427 fData->Divide(h2,h1);
431 //____________________________________________________________________
432 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2, Option_t *option)
435 //Divide aGrid1 and aGrid2 and deposit the result into the current one
436 //binomial errors are supported
439 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
440 AliError("Different number of variables, cannot divide the grids");
444 if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
446 THnSparse *h1= aGrid1->GetGrid();
447 THnSparse *h2= aGrid2->GetGrid();
448 fData->Divide(h1,h2,c1,c2,option);
452 //____________________________________________________________________
453 void AliCFGridSparse::Rebin(const Int_t* group)
456 // rebin the grid according to Rebin() as in THnSparse
457 // Please notice that the original number of bins on
458 // a given axis has to be divisible by the rebin group.
461 for(Int_t i=0;i<GetNVar();i++){
462 if (group[i]!=1) AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
465 THnSparse *rebinned =fData->Rebin(group);
469 //____________________________________________________________________
470 void AliCFGridSparse::Scale(Long_t index, const Double_t *fact)
473 //scale content of a certain cell by (positive) fact (with error)
476 if (GetElement(index)==0 || fact[0]==0) return;
478 Double_t in[2], out[2];
479 in[0]=GetElement(index);
480 in[1]=GetElementError(index);
481 GetScaledValues(fact,in,out);
482 SetElement(index,out[0]);
483 if (fSumW2) SetElementError(index,out[1]);
485 //____________________________________________________________________
486 void AliCFGridSparse::Scale(const Int_t *bin, const Double_t *fact)
489 //scale content of a certain cell by (positive) fact (with error)
491 if(GetElement(bin)==0 || fact[0]==0)return;
493 Double_t in[2], out[2];
494 in[0]=GetElement(bin);
495 in[1]=GetElementError(bin);
496 GetScaledValues(fact,in,out);
497 SetElement(bin,out[0]);
498 if(fSumW2)SetElementError(bin,out[1]);
501 //____________________________________________________________________
502 void AliCFGridSparse::Scale(const Double_t *var, const Double_t *fact)
505 //scale content of a certain cell by (positive) fact (with error)
507 if(GetElement(var)==0 || fact[0]==0)return;
509 Double_t in[2], out[2];
510 in[0]=GetElement(var);
511 in[1]=GetElementError(var);
512 GetScaledValues(fact,in,out);
513 SetElement(var,out[0]);
514 if(fSumW2)SetElementError(var,out[1]);
517 //____________________________________________________________________
518 void AliCFGridSparse::Scale(const Double_t* fact)
521 //scale contents of the whole grid by fact
524 for (Long_t iel=0; iel<GetNFilledBins(); iel++) {
528 //____________________________________________________________________
529 Long_t AliCFGridSparse::GetEmptyBins() const {
534 return (GetNBinsTotal() - GetNFilledBins()) ;
537 //____________________________________________________________________
538 Int_t AliCFGridSparse::CheckStats(Double_t thr) const
541 // Count the cells below a certain threshold
544 for (Int_t i=0; i<GetNBinsTotal(); i++) {
545 if (GetElement(i)<thr) ncellsLow++;
550 //_____________________________________________________________________
551 Double_t AliCFGridSparse::GetIntegral() const
556 return fData->ComputeIntegral();
559 //____________________________________________________________________
560 Long64_t AliCFGridSparse::Merge(TCollection* list)
563 // Merge a list of AliCFGridSparse with this (needed for PROOF).
564 // Returns the number of merged objects (including this).
573 TIterator* iter = list->MakeIterator();
577 while ((obj = iter->Next())) {
578 AliCFGridSparse* entry = dynamic_cast<AliCFGridSparse*> (obj);
588 //____________________________________________________________________
589 void AliCFGridSparse::GetScaledValues(const Double_t *fact, const Double_t *in, Double_t *out) const{
591 // scale input *in and its error by (positive) fact (with error)
592 // and erite it to *out
594 out[0]=in[0]*fact[0];
595 out[1]=TMath::Sqrt(in[1]*in[1]/in[0]/in[0]
596 +fact[1]*fact[1]/fact[0]/fact[0])*out[0];
600 //____________________________________________________________________
601 void AliCFGridSparse::Copy(TObject& c) const
607 AliCFGridSparse& target = (AliCFGridSparse &) c;
608 target.fSumW2 = fSumW2 ;
610 target.fData = (THnSparse*)fData->Clone();
614 //____________________________________________________________________
615 TH1* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
618 // return a slice on variables iVar1 (and optionnally iVar2 (and iVar3)) while axis ranges are defined with varMin,varMax
619 // arrays varMin and varMax contain the min and max values of each variable.
620 // therefore varMin and varMax must have their dimensions equal to GetNVar()
621 // If useBins=true, varMin and varMax are taken as bin numbers
622 // if varmin or varmax point to null, all the range is taken, including over- and underflows
624 THnSparse* clone = (THnSparse*)fData->Clone();
625 if (varMin != 0x0 && varMax != 0x0) {
626 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) SetAxisRange(clone->GetAxis(iAxis),varMin[iAxis],varMax[iAxis],useBins);
629 TH1* projection = 0x0 ;
631 GetProjectionName (name ,iVar1,iVar2,iVar3);
632 GetProjectionTitle(title,iVar1,iVar2,iVar3);
636 if (iVar1 >= GetNVar() || iVar1 < 0 ) {
637 AliError("Non-existent variable, return NULL");
640 projection = (TH1D*)clone->Projection(iVar1);
641 projection->SetTitle(Form("%s_proj-%s",GetTitle(),GetVarTitle(iVar1)));
642 for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
643 Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
644 TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
645 if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
649 if (iVar1 >= GetNVar() || iVar1 < 0 ||
650 iVar2 >= GetNVar() || iVar2 < 0 ) {
651 AliError("Non-existent variable, return NULL");
654 projection = (TH2D*)clone->Projection(iVar2,iVar1);
655 for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
656 Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
657 TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
658 if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
660 for (Int_t iBin=1; iBin<=projection->GetNbinsY(); iBin++) {
661 Int_t origBin = GetAxis(iVar2)->GetFirst()+iBin-1;
662 TString binLabel = GetAxis(iVar2)->GetBinLabel(origBin) ;
663 if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
668 if (iVar1 >= GetNVar() || iVar1 < 0 ||
669 iVar2 >= GetNVar() || iVar2 < 0 ||
670 iVar3 >= GetNVar() || iVar3 < 0 ) {
671 AliError("Non-existent variable, return NULL");
674 projection = (TH3D*)clone->Projection(iVar1,iVar2,iVar3);
675 for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
676 Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
677 TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
678 if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
680 for (Int_t iBin=1; iBin<=projection->GetNbinsY(); iBin++) {
681 Int_t origBin = GetAxis(iVar2)->GetFirst()+iBin-1;
682 TString binLabel = GetAxis(iVar2)->GetBinLabel(origBin) ;
683 if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
685 for (Int_t iBin=1; iBin<=projection->GetNbinsZ(); iBin++) {
686 Int_t origBin = GetAxis(iVar3)->GetFirst()+iBin-1;
687 TString binLabel = GetAxis(iVar3)->GetBinLabel(origBin) ;
688 if (binLabel.CompareTo("") != 0) projection->GetZaxis()->SetBinLabel(iBin,binLabel);
692 projection->SetName (name .Data());
693 projection->SetTitle(title.Data());
699 //____________________________________________________________________
700 void AliCFGridSparse::SetAxisRange(TAxis* axis, Double_t min, Double_t max, Bool_t useBins) const {
702 // sets axis range, and forces bit TAxis::kAxisRange
705 if (useBins) axis->SetRange ((Int_t)min,(Int_t)max);
706 else axis->SetRangeUser( min, max);
707 //axis->SetBit(TAxis::kAxisRange); // uncomment when ROOT TAxis is fixed
710 //____________________________________________________________________
711 void AliCFGridSparse::SetRangeUser(Int_t iVar, Double_t varMin, Double_t varMax, Bool_t useBins) const {
713 // set range of axis iVar.
715 SetAxisRange(fData->GetAxis(iVar),varMin,varMax,useBins);
716 //AliInfo(Form("AliCFGridSparse axis %d range has been modified",iVar));
717 TAxis* currAxis = fData->GetAxis(iVar);
718 TString outString = Form("%s new range: %.1f < %s < %.1f", GetName(), currAxis->GetBinLowEdge(currAxis->GetFirst()), currAxis->GetTitle(), currAxis->GetBinUpEdge(currAxis->GetLast()));
719 TString binLabel = currAxis->GetBinLabel(currAxis->GetFirst());
720 if ( ! binLabel.IsNull() ) {
722 for ( Int_t ibin = currAxis->GetFirst(); ibin <= currAxis->GetLast(); ibin++ ) {
723 outString += Form("%s ", currAxis->GetBinLabel(ibin));
727 AliWarning(outString.Data());
731 //____________________________________________________________________
732 void AliCFGridSparse::SetRangeUser(const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const {
734 // set range of every axis. varMin and varMax must be of dimension GetNVar()
736 for (Int_t iAxis=0; iAxis<GetNVar() ; iAxis++) { // set new range for every axis
737 SetRangeUser(iAxis,varMin[iAxis],varMax[iAxis], useBins);
739 AliInfo("AliCFGridSparse axes ranges have been modified");
742 //____________________________________________________________________
743 Float_t AliCFGridSparse::GetOverFlows(Int_t ivar, Bool_t exclusive) const
746 // Returns overflows in variable ivar
747 // Set 'exclusive' to true for an exclusive check on variable ivar
749 Int_t* bin = new Int_t[GetNVar()];
750 memset(bin, 0, sizeof(Int_t) * GetNVar());
752 for (Long64_t i = 0; i < fData->GetNbins(); i++) {
753 Double_t v = fData->GetBinContent(i, bin);
756 for(Int_t j=0;j<GetNVar();j++){
758 if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
761 if(bin[ivar]==GetNBins(ivar)+1 && add) ovfl+=v;
768 //____________________________________________________________________
769 Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar, Bool_t exclusive) const
772 // Returns exclusive overflows in variable ivar
773 // Set 'exclusive' to true for an exclusive check on variable ivar
775 Int_t* bin = new Int_t[GetNVar()];
776 memset(bin, 0, sizeof(Int_t) * GetNVar());
778 for (Long64_t i = 0; i < fData->GetNbins(); i++) {
779 Double_t v = fData->GetBinContent(i, bin);
782 for(Int_t j=0;j<GetNVar();j++){
784 if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
787 if(bin[ivar]==0 && add) unfl+=v;
795 //____________________________________________________________________
796 void AliCFGridSparse::Smooth() {
798 // smoothing function: TO USE WITH CARE
801 AliInfo("Your GridSparse is going to be smoothed");
802 AliInfo(Form("N TOTAL BINS : %li",GetNBinsTotal()));
803 AliInfo(Form("N FILLED BINS : %li",GetNFilledBins()));
804 AliCFUnfolding::SmoothUsingNeighbours(fData);