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