Patch from Renaud for the correction framework
[u/mrichter/AliRoot.git] / CORRFW / AliCFGridSparse.cxx
CommitLineData
db6722a5 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"
7036630f 36#include "AliCFUnfolding.h"
db6722a5 37
38//____________________________________________________________________
39ClassImp(AliCFGridSparse)
40
41//____________________________________________________________________
42AliCFGridSparse::AliCFGridSparse() :
fb494025 43 AliCFFrame(),
44 fSumW2(kFALSE),
db6722a5 45 fData(0x0)
46{
47 // default constructor
48}
49//____________________________________________________________________
50AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
fb494025 51 AliCFFrame(name,title),
52 fSumW2(kFALSE),
db6722a5 53 fData(0x0)
54{
55 // default constructor
56}
57//____________________________________________________________________
fb494025 58AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, Int_t nVarIn, const Int_t * nBinIn) :
59 AliCFFrame(name,title),
60 fSumW2(kFALSE),
db6722a5 61 fData(0x0)
62{
63 //
64 // main constructor
65 //
66
fb494025 67 fData = new THnSparseF(name,title,nVarIn,nBinIn);
db6722a5 68}
fb494025 69
db6722a5 70//____________________________________________________________________
fb494025 71AliCFGridSparse::~AliCFGridSparse()
db6722a5 72{
73 //
fb494025 74 // destructor
db6722a5 75 //
fb494025 76 if (fData) delete fData;
db6722a5 77}
78
79//____________________________________________________________________
fb494025 80AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
81 AliCFFrame(c),
82 fSumW2(kFALSE),
83 fData(0x0)
db6722a5 84{
85 //
fb494025 86 // copy constructor
db6722a5 87 //
fb494025 88 ((AliCFGridSparse &)c).Copy(*this);
db6722a5 89}
90
91//____________________________________________________________________
fb494025 92AliCFGridSparse& AliCFGridSparse::operator=(const AliCFGridSparse &c)
db6722a5 93{
94 //
95 // assigment operator
96 //
fb494025 97 if (this != &c) c.Copy(*this);
db6722a5 98 return *this;
99}
100
101//____________________________________________________________________
aa29aca9 102void 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//____________________________________________________________________
fb494025 115void AliCFGridSparse::SetBinLimits(Int_t ivar, const Double_t *array)
db6722a5 116{
117 //
118 // setting the arrays containing the bin limits
119 //
120 fData->SetBinEdges(ivar, array);
db6722a5 121}
122
123//____________________________________________________________________
fb494025 124void AliCFGridSparse::Fill(const Double_t *var, Double_t weight)
db6722a5 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//___________________________________________________________________
135TH1D *AliCFGridSparse::Project(Int_t ivar) const
136{
137 //
138 // Make a 1D projection along variable ivar
139 //
140
db6722a5 141 TH1D *hist=fData->Projection(ivar);
fb494025 142 hist->SetXTitle(GetVarTitle(ivar));
91e32eeb 143 hist->SetName(Form("%s_proj-%s",GetName(),GetVarTitle(ivar)));
fb494025 144 hist->SetTitle(Form("%s: projection on %s",GetTitle(),GetVarTitle(ivar)));
b12ab491 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 }
db6722a5 149 return hist;
db6722a5 150}
151//___________________________________________________________________
152TH2D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2) const
153{
154 //
155 // Make a 2D projection along variables ivar1 & ivar2
156 //
157
318f64b1 158 TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
fb494025 159 hist->SetXTitle(GetVarTitle(ivar1));
160 hist->SetYTitle(GetVarTitle(ivar2));
91e32eeb 161 hist->SetName(Form("%s_proj-%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
fb494025 162 hist->SetTitle(Form("%s: projection on %s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
b12ab491 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 }
db6722a5 171 return hist;
db6722a5 172}
173//___________________________________________________________________
174TH3D *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 //
db6722a5 179
180 TH3D *hist=fData->Projection(ivar1,ivar2,ivar3);
fb494025 181 hist->SetXTitle(GetVarTitle(ivar1));
182 hist->SetYTitle(GetVarTitle(ivar2));
183 hist->SetZTitle(GetVarTitle(ivar3));
91e32eeb 184 hist->SetName(Form("%s_proj-%s,%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
fb494025 185 hist->SetTitle(Form("%s: projection on %s-%s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
b12ab491 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 }
db6722a5 198 return hist;
db6722a5 199}
200
9105291d 201//___________________________________________________________________
98a5f772 202AliCFGridSparse* AliCFGridSparse::Project(Int_t nVars, const Int_t* vars, const Double_t* varMin, const Double_t* varMax, Bool_t useBins) const
9105291d 203{
38b1447f 204 //
205 // projects the grid on the nVars dimensions defined in vars.
206 // axis ranges can be defined in arrays varMin, varMax
98a5f772 207 // If useBins=true, varMin and varMax are taken as bin numbers
38b1447f 208 //
9105291d 209
210 // binning for new grid
211 Int_t* bins = new Int_t[nVars];
212 for (Int_t iVar=0; iVar<nVars; iVar++) {
fb494025 213 bins[iVar] = GetNBins(vars[iVar]);
9105291d 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());
38b1447f 221 if (varMin && varMax) {
fb494025 222 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
98a5f772 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]);
38b1447f 225 }
9105291d 226 }
38b1447f 227 else AliInfo("Keeping same axis ranges");
228
9105291d 229 out->SetGrid(clone->Projection(nVars,vars));
4b7819c4 230 delete [] bins;
9105291d 231 return out;
232}
233
fb494025 234
db6722a5 235//____________________________________________________________________
fb494025 236Float_t AliCFGridSparse::GetBinCenter(Int_t ivar, Int_t ibin) const
db6722a5 237{
238 //
fb494025 239 // Returns the center of specified bin for variable axis ivar
240 //
241
242 return (Float_t) fData->GetAxis(ivar)->GetBinCenter(ibin);
db6722a5 243}
244
245//____________________________________________________________________
fb494025 246Float_t AliCFGridSparse::GetBinSize(Int_t ivar, Int_t ibin) const
db6722a5 247{
248 //
fb494025 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);
db6722a5 253}
254
db6722a5 255//____________________________________________________________________
256Float_t AliCFGridSparse::GetEntries() const
257{
258 //
259 // total entries (including overflows and underflows)
260 //
261
262 return fData->GetEntries();
263}
264
265//____________________________________________________________________
fb494025 266Float_t AliCFGridSparse::GetElement(Long_t index) const
db6722a5 267{
268 //
fb494025 269 // Returns content of grid element index
db6722a5 270 //
db6722a5 271
fb494025 272 return fData->GetBinContent(index);
db6722a5 273}
274//____________________________________________________________________
fb494025 275Float_t AliCFGridSparse::GetElement(const Int_t *bin) const
db6722a5 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//____________________________________________________________________
fb494025 284Float_t AliCFGridSparse::GetElement(const Double_t *var) const
db6722a5 285{
286 //
287 // Get the content in a bin corresponding to a set of input variables
288 //
289
fb494025 290 Long_t index = fData->GetBin(var,kFALSE);
291 if (index<0) return 0.;
292 return fData->GetBinContent(index);
db6722a5 293}
294
295//____________________________________________________________________
fb494025 296Float_t AliCFGridSparse::GetElementError(Long_t index) const
db6722a5 297{
298 //
fb494025 299 // Returns the error on the content
db6722a5 300 //
301
fb494025 302 return fData->GetBinContent(index);
db6722a5 303}
304//____________________________________________________________________
fb494025 305Float_t AliCFGridSparse::GetElementError(const Int_t *bin) const
db6722a5 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//____________________________________________________________________
fb494025 314Float_t AliCFGridSparse::GetElementError(const Double_t *var) const
db6722a5 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)
fb494025 321 if (index<0) return 0.;
322 return fData->GetBinError(index);
db6722a5 323}
324
db6722a5 325//____________________________________________________________________
fb494025 326void AliCFGridSparse::SetElement(Long_t index, Float_t val)
db6722a5 327{
328 //
fb494025 329 // Sets grid element value
db6722a5 330 //
fb494025 331 Int_t* bin = new Int_t[GetNVar()];
332 fData->GetBinContent(index,bin); //affects the bin coordinates
333 SetElement(bin,val);
334 delete [] bin ;
db6722a5 335}
fb494025 336
db6722a5 337//____________________________________________________________________
fb494025 338void AliCFGridSparse::SetElement(const Int_t *bin, Float_t val)
db6722a5 339{
340 //
341 // Sets grid element of bin indeces bin to val
342 //
343 fData->SetBinContent(bin,val);
344}
345//____________________________________________________________________
fb494025 346void AliCFGridSparse::SetElement(const Double_t *var, Float_t val)
db6722a5 347{
348 //
349 // Set the content in a bin to value val corresponding to a set of input variables
350 //
fb494025 351 Long_t index=fData->GetBin(var,kTRUE); //THnSparse index: allocate the cell
352 Int_t *bin = new Int_t[GetNVar()];
db6722a5 353 fData->GetBinContent(index,bin); //trick to access the array of bins
fb494025 354 SetElement(bin,val);
db6722a5 355 delete [] bin;
db6722a5 356}
357
358//____________________________________________________________________
fb494025 359void AliCFGridSparse::SetElementError(Long_t index, Float_t val)
db6722a5 360{
361 //
362 // Sets grid element iel error to val (linear indexing) in AliCFFrame
363 //
fb494025 364 Int_t *bin = new Int_t[GetNVar()];
365 fData->GetBinContent(index,bin);
366 SetElementError(bin,val);
db6722a5 367 delete [] bin;
368}
fb494025 369
db6722a5 370//____________________________________________________________________
fb494025 371void AliCFGridSparse::SetElementError(const Int_t *bin, Float_t val)
db6722a5 372{
373 //
374 // Sets grid element error of bin indeces bin to val
375 //
376 fData->SetBinError(bin,val);
377}
378//____________________________________________________________________
fb494025 379void AliCFGridSparse::SetElementError(const Double_t *var, Float_t val)
db6722a5 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
fb494025 385 Int_t *bin = new Int_t[GetNVar()];
db6722a5 386 fData->GetBinContent(index,bin); //trick to access the array of bins
fb494025 387 SetElementError(bin,val);
db6722a5 388 delete [] bin;
389}
390
391//____________________________________________________________________
392void AliCFGridSparse::SumW2()
393{
394 //
395 //set calculation of the squared sum of the weighted entries
396 //
397 if(!fSumW2){
398 fData->CalculateErrors(kTRUE);
399 }
db6722a5 400 fSumW2=kTRUE;
401}
402
403//____________________________________________________________________
fb494025 404void AliCFGridSparse::Add(const AliCFGridSparse* aGrid, Double_t c)
db6722a5 405{
406 //
407 //add aGrid to the current one
408 //
409
fb494025 410 if (aGrid->GetNVar() != GetNVar()){
411 AliError("Different number of variables, cannot add the grids");
db6722a5 412 return;
413 }
db6722a5 414
fb494025 415 if (!fSumW2 && aGrid->GetSumW2()) SumW2();
416 fData->Add(aGrid->GetGrid(),c);
db6722a5 417}
418
419//____________________________________________________________________
fb494025 420void AliCFGridSparse::Add(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
db6722a5 421{
422 //
423 //Add aGrid1 and aGrid2 and deposit the result into the current one
424 //
425
fb494025 426 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
db6722a5 427 AliInfo("Different number of variables, cannot add the grids");
428 return;
429 }
db6722a5 430
fb494025 431 if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
db6722a5 432
433 fData->Reset();
fb494025 434 fData->Add(aGrid1->GetGrid(),c1);
435 fData->Add(aGrid2->GetGrid(),c2);
db6722a5 436}
437
438//____________________________________________________________________
fb494025 439void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid, Double_t c)
db6722a5 440{
441 //
442 // Multiply aGrid to the current one
443 //
444
fb494025 445 if (aGrid->GetNVar() != GetNVar()) {
446 AliError("Different number of variables, cannot multiply the grids");
db6722a5 447 return;
448 }
db6722a5 449
fb494025 450 if(!fSumW2 && aGrid->GetSumW2()) SumW2();
451 THnSparse *h = aGrid->GetGrid();
db6722a5 452 fData->Multiply(h);
453 fData->Scale(c);
db6722a5 454}
455
456//____________________________________________________________________
fb494025 457void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
db6722a5 458{
459 //
460 //Multiply aGrid1 and aGrid2 and deposit the result into the current one
461 //
462
fb494025 463 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
464 AliError("Different number of variables, cannot multiply the grids");
db6722a5 465 return;
fb494025 466 }
db6722a5 467
fb494025 468 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
db6722a5 469
470 fData->Reset();
fb494025 471 THnSparse *h1 = aGrid1->GetGrid();
472 THnSparse *h2 = aGrid2->GetGrid();
db6722a5 473 h2->Multiply(h1);
474 h2->Scale(c1*c2);
475 fData->Add(h2);
476}
477
db6722a5 478//____________________________________________________________________
fb494025 479void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid, Double_t c)
db6722a5 480{
481 //
482 // Divide aGrid to the current one
483 //
484
fb494025 485 if (aGrid->GetNVar() != GetNVar()) {
486 AliError("Different number of variables, cannot divide the grids");
db6722a5 487 return;
488 }
db6722a5 489
fb494025 490 if (!fSumW2 && aGrid->GetSumW2()) SumW2();
db6722a5 491
d7c1cb54 492 THnSparse *h1 = aGrid->GetGrid();
493 THnSparse *h2 = (THnSparse*)fData->Clone();
494 fData->Divide(h2,h1);
db6722a5 495 fData->Scale(c);
db6722a5 496}
497
498//____________________________________________________________________
fb494025 499void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2, Option_t *option)
db6722a5 500{
501 //
502 //Divide aGrid1 and aGrid2 and deposit the result into the current one
fb494025 503 //binomial errors are supported
db6722a5 504 //
505
fb494025 506 if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
507 AliError("Different number of variables, cannot divide the grids");
db6722a5 508 return;
509 }
db6722a5 510
fb494025 511 if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
db6722a5 512
fb494025 513 THnSparse *h1= aGrid1->GetGrid();
514 THnSparse *h2= aGrid2->GetGrid();
db6722a5 515 fData->Divide(h1,h2,c1,c2,option);
516}
517
518
519//____________________________________________________________________
7411edfd 520void 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
fb494025 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]));
7411edfd 530 }
531
532 THnSparse *rebinned =fData->Rebin(group);
533 fData->Reset();
534 fData = rebinned;
fb494025 535}
536//____________________________________________________________________
537void AliCFGridSparse::Scale(Long_t index, const Double_t *fact)
538{
539 //
540 //scale content of a certain cell by (positive) fact (with error)
541 //
7411edfd 542
fb494025 543 if (GetElement(index)==0 || fact[0]==0) return;
7411edfd 544
fb494025 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//____________________________________________________________________
553void 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;
7411edfd 559
fb494025 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//____________________________________________________________________
569void 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;
7411edfd 575
fb494025 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//____________________________________________________________________
585void 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);
7411edfd 593 }
fb494025 594}
595//____________________________________________________________________
596Long_t AliCFGridSparse::GetEmptyBins() const {
597 //
598 // Get empty bins
599 //
7411edfd 600
fb494025 601 return (GetNBinsTotal() - GetNFilledBins()) ;
602}
7411edfd 603
fb494025 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//____________________________________________________________________
646Int_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}
7411edfd 657
fb494025 658//_____________________________________________________________________
659Double_t AliCFGridSparse::GetIntegral() const
660{
661 //
662 // Get full Integral
663 //
664 return fData->ComputeIntegral();
665}
7411edfd 666
fb494025 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//____________________________________________________________________
753Long64_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;
7411edfd 762
fb494025 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;
7411edfd 779}
fb494025 780
781//____________________________________________________________________
782void 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
7411edfd 793//____________________________________________________________________
db6722a5 794void AliCFGridSparse::Copy(TObject& c) const
795{
796 //
797 // copy function
798 //
fb494025 799 AliCFFrame::Copy(c);
db6722a5 800 AliCFGridSparse& target = (AliCFGridSparse &) c;
fb494025 801 target.fSumW2 = fSumW2 ;
802 if (fData) {
803 target.fData = (THnSparse*)fData->Clone();
804 }
db6722a5 805}
806
c8df672e 807//____________________________________________________________________
98a5f772 808TH1D* AliCFGridSparse::Slice(Int_t iVar, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
c8df672e 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.
fb494025 813 // therefore varMin and varMax must have their dimensions equal to GetNVar()
98a5f772 814 // If useBins=true, varMin and varMax are taken as bin numbers
c8df672e 815
816 THnSparse* clone = (THnSparse*)fData->Clone();
fb494025 817 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
98a5f772 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]);
c8df672e 820 }
b12ab491 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 }
53ab59fe 836 return projection ;
c8df672e 837}
838
839//____________________________________________________________________
98a5f772 840TH2D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
c8df672e 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.
fb494025 845 // therefore varMin and varMax must have their dimensions equal to GetNVar()
98a5f772 846 // If useBins=true, varMin and varMax are taken as bin numbers
c8df672e 847
848 THnSparse* clone = (THnSparse*)fData->Clone();
fb494025 849 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
98a5f772 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]);
c8df672e 852 }
b12ab491 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 }
53ab59fe 872 return projection ;
c8df672e 873}
874
875//____________________________________________________________________
98a5f772 876TH3D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
c8df672e 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.
fb494025 881 // therefore varMin and varMax must have their dimensions equal to GetNVar()
98a5f772 882 // If useBins=true, varMin and varMax are taken as bin numbers
c8df672e 883
884 THnSparse* clone = (THnSparse*)fData->Clone();
fb494025 885 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
98a5f772 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]);
c8df672e 888 }
b12ab491 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 }
53ab59fe 912 return projection ;
c8df672e 913}
914
915//____________________________________________________________________
9105291d 916void 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//____________________________________________________________________
fb494025 925void AliCFGridSparse::SetRangeUser(const Double_t *varMin, const Double_t *varMax) {
c8df672e 926 //
fb494025 927 // set range of every axis. varMin and varMax must be of dimension GetNVar()
c8df672e 928 //
fb494025 929 for (Int_t iAxis=0; iAxis<GetNVar() ; iAxis++) { // set new range for every axis
9105291d 930 SetRangeUser(iAxis,varMin[iAxis],varMax[iAxis]);
c8df672e 931 }
9105291d 932 AliWarning("THnSparse axes ranges have been modified");
933}
934
935//____________________________________________________________________
936void AliCFGridSparse::UseAxisRange(Bool_t b) const {
fb494025 937 for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) fData->GetAxis(iAxis)->SetBit(TAxis::kAxisRange,b);
938}
939
940//____________________________________________________________________
941Float_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//____________________________________________________________________
967Float_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;
c8df672e 990}
fb494025 991
7036630f 992
993//____________________________________________________________________
994void 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}