]> git.uio.no Git - u/mrichter/AliRoot.git/blame - CORRFW/AliCFGrid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update Container Classes to be able to handle also THnSparse-based implementation...
[u/mrichter/AliRoot.git] / CORRFW / AliCFGrid.cxx
CommitLineData
563113d0 1/* $Id$ */
1e9dad92 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16//---------------------------------------------------------------------//
17// //
18// AliCFGrid Class //
19// Class to accumulate data on an N-dimensional grid, to be used //
20// as input to get corrections for Reconstruction & Trigger efficiency //
21// The class uses a one-dimensional array of floats to store the grid //
22// --Author : S.Arcelli //
23// Still to be done: //
24// --Implement methods to merge cells //
25// --Interpolate among bins in a range //
26// This implementation will be aventually replaced byAliCFGridSparse //
27//---------------------------------------------------------------------//
563113d0 28//
29//
30#include <AliLog.h>
31#include "AliCFGrid.h"
32#include "TMath.h"
33#include "TROOT.h"
1e9dad92 34#include "TH1D.h"
35#include "TH2D.h"
36#include "TH3D.h"
563113d0 37
38//____________________________________________________________________
39ClassImp(AliCFGrid)
40
41//____________________________________________________________________
42AliCFGrid::AliCFGrid() :
1e9dad92 43 AliCFVGrid(),
563113d0 44 fNunflTot(0),
45 fNovflTot(0),
46 fNentriesTot(0),
47 fNunfl(0x0),
48 fNovfl(0x0),
49 fData(0x0),
50 fErr2(0x0)
51{
52 // default constructor
53}
54//____________________________________________________________________
55AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title) :
1e9dad92 56 AliCFVGrid(name,title),
563113d0 57 fNunflTot(0),
58 fNovflTot(0),
59 fNentriesTot(0),
60 fNunfl(0x0),
61 fNovfl(0x0),
62 fData(0x0),
63 fErr2(0x0)
64{
65 // default constructor
66}
67
68//____________________________________________________________________
1e9dad92 69AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title, const Int_t nVarIn, const Int_t * nBinIn, const Double_t *binLimitsIn) :
70 AliCFVGrid(name,title,nVarIn,nBinIn,binLimitsIn),
563113d0 71 fNunflTot(0),
72 fNovflTot(0),
73 fNentriesTot(0),
74 fNunfl(0x0),
75 fNovfl(0x0),
76 fData(0x0),
77 fErr2(0x0)
78{
79 //
80 // main constructor
81 //
1e9dad92 82
83 //The over/underflows
563113d0 84 fNunfl=new Float_t[fNVar];
85 fNovfl= new Float_t[fNVar];
86
87
563113d0 88 //Initialization
563113d0 89 fNunflTot =0;
90 fNovflTot =0;
91 fNentriesTot =0;
92 for(Int_t j=0;j<fNVar;j++){
93 fNunfl[j] =0;
94 fNovfl[j] =0;
95 }
1e9dad92 96
97
98 // the grid
99
100 fData = new Float_t[fNDim]; //num
101
102 //Initialization
103
563113d0 104}
105
106//____________________________________________________________________
107AliCFGrid::AliCFGrid(const AliCFGrid& c) :
1e9dad92 108 AliCFVGrid(c),
563113d0 109 fNunflTot(0),
110 fNovflTot(0),
111 fNentriesTot(0),
112 fNunfl(0x0),
113 fNovfl(0x0),
114 fData(0x0),
115 fErr2(0x0)
116{
117 //
118 // copy constructor
119 //
120 ((AliCFGrid &)c).Copy(*this);
121}
122
123//____________________________________________________________________
124AliCFGrid::~AliCFGrid()
125{
126 //
127 // destructor
128 //
1e9dad92 129 if(fNunfl)delete fNunfl;
130 if(fNovfl)delete fNovfl;
131 if(fData)delete fData;
563113d0 132 if(fSumW2)delete fErr2;
563113d0 133
134}
563113d0 135//____________________________________________________________________
136AliCFGrid &AliCFGrid::operator=(const AliCFGrid &c)
137{
138 //
139 // assigment operator
140 //
141 if (this != &c)
142 ((AliCFGrid &) c).Copy(*this);
143 return *this;
144}
145//____________________________________________________________________
146Float_t AliCFGrid::GetElement(Int_t bin) const
147{
148 //
149 // Returns grid element bin
150 //
151 if(bin>=fNDim){
152 AliInfo(Form(" element index outside the grid, return -1"));
153 return -1.;
154 }
155 return fData[bin];
156}
157//____________________________________________________________________
158Float_t AliCFGrid::GetElement(Int_t *bin) const
159{
160 //
161 // Get the content in a bin corresponding to a set of bin indexes
162 //
1e9dad92 163 //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
164 for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
165 Int_t ind =GetBinIndex(fIndex);
166 return GetElement(ind);
563113d0 167}
168//____________________________________________________________________
1e9dad92 169Float_t AliCFGrid::GetElement(Double_t *var) const
563113d0 170{
171 //
172 // Get the content in a bin corresponding to a set of input variables
173 //
174 Int_t unfl=0;
175 Int_t ovfl=0;
176 for(Int_t i=0;i<fNVar;i++){
177 Int_t nbins=fNVarBins[i]+1;
178 Float_t *bins=new Float_t[nbins];
179 for(Int_t ibin =0;ibin<nbins;ibin++){
180 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
181 }
182
183 fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
184
185 //underflows
186
187 if(var[i] < bins[0]){
188 unfl=1;
189 }
190
191 //overflows
192
193 if(var[i] > bins[nbins-1]){
194 ovfl=1;
195 }
196 delete [] bins;
197 }
198
1e9dad92 199
200 //move to the TH/THnSparse convention in N-dim bin numbering
201 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
202
563113d0 203 if(!(ovfl==1 || unfl==1)){
204 return GetElement(fIndex);
205 }
206 else{
207 AliInfo(Form(" input variables outside the grid, return -1"));
208 return -1.;
209 }
210}
211//____________________________________________________________________
212Float_t AliCFGrid::GetElementError(Int_t iel) const
213{
214 //
1e9dad92 215 // Return the error on grid element iel
563113d0 216 //
217 if(iel>=fNDim){
218 AliInfo(Form(" element index outside the grid, return -1"));
219 return -1.;
220 }
1e9dad92 221 if(fSumW2)return TMath::Sqrt(fErr2[iel]);
222 return TMath::Sqrt(fData[iel]);
563113d0 223}
224//____________________________________________________________________
225Float_t AliCFGrid::GetElementError(Int_t *bin) const
226{
227 //
1e9dad92 228 // Get the error in a bin corresponding to a set of bin indeces
563113d0 229 //
1e9dad92 230 //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
231 for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
232 Int_t ind =GetBinIndex(fIndex);
233 return GetElementError(ind);
563113d0 234
235}
236//____________________________________________________________________
1e9dad92 237Float_t AliCFGrid::GetElementError(Double_t *var) const
563113d0 238{
239 //
1e9dad92 240 // Get the error in a bin corresponding to a set of input variables
563113d0 241 //
242 Int_t unfl=0;
243 Int_t ovfl=0;
244 for(Int_t i=0;i<fNVar;i++){
245 Int_t nbins=fNVarBins[i]+1;
246 Float_t *bins=new Float_t[nbins];
247 for(Int_t ibin =0;ibin<nbins;ibin++){
248 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
249 }
250
251 fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
252 //underflows
253
254 if(var[i] < bins[0]){
255 unfl=1;
256 }
257
258 //overflows
259
260 if(var[i] > bins[nbins-1]){
261 ovfl=1;
262 }
263 delete [] bins;
264 }
265
1e9dad92 266 //move to the TH/THnSparse convention in N-dim bin numbering
267 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
563113d0 268 if(!(ovfl==1 || unfl==1)){
269 return GetElementError(fIndex);
270 }
271 else{
272 AliInfo(Form(" input variables outside the grid, return -1"));
273 return -1.;
274 }
275}
276//____________________________________________________________________
277void AliCFGrid::SetElement(Int_t iel, Float_t val)
278{
279 //
280 // Sets grid element iel to val
281 //
282 if(iel>=fNDim){
283 AliInfo(Form(" element index outside the grid, no value set"));
284 }else {
285 fData[iel]=val;
286 }
287}
288//____________________________________________________________________
289void AliCFGrid::SetElement(Int_t *bin, Float_t val)
290{
291 //
292 // Sets grid element of bin indeces bin to val
293 //
1e9dad92 294 //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
295 for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
296 Int_t ind =GetBinIndex(fIndex);
297 SetElement(ind,val);
563113d0 298}
299//____________________________________________________________________
1e9dad92 300void AliCFGrid::SetElement(Double_t *var, Float_t val)
563113d0 301{
302 //
303 // Set the content in a bin to value val corresponding to a set of input variables
304 //
305 Int_t unfl=0;
306 Int_t ovfl=0;
307 for(Int_t i=0;i<fNVar;i++){
308 Int_t nbins=fNVarBins[i]+1;
309 Float_t *bins=new Float_t[nbins];
310 for(Int_t ibin =0;ibin<nbins;ibin++){
311 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
312 }
313
314 fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
315 //underflows
316
317 if(var[i] < bins[0]){
318 unfl=1;
319 }
320
321 //overflows
322
323 if(var[i] > bins[nbins-1]){
324 ovfl=1;
325 }
326 delete [] bins;
327 }
328
1e9dad92 329 //move to the TH/THnSparse convention in N-dim bin numbering
330 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
563113d0 331 if(!(ovfl==1 || unfl==1)){
332 SetElement(fIndex,val);
333 }
334 else{
335 AliInfo(Form(" input variables outside the grid, no value set"));
336 }
337}
338//____________________________________________________________________
339void AliCFGrid::SetElementError(Int_t iel, Float_t val)
340{
341 //
1e9dad92 342 // Set squared error on grid element iel to val*val
563113d0 343 //
344 if(iel>=fNDim){
345 AliInfo(Form(" element index outside the grid, no value set"));
346 return;
347 }
348 if(!fErr2)SumW2();
1e9dad92 349 fErr2[iel]=val*val;
563113d0 350}
351//____________________________________________________________________
352void AliCFGrid::SetElementError(Int_t *bin, Float_t val)
353{
354 //
355 // Set squared error to val on grid element of bin indeces bin
356 //
1e9dad92 357 //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
358 for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
359 Int_t ind =GetBinIndex(fIndex);
360 SetElementError(ind,val);
563113d0 361}
362//____________________________________________________________________
1e9dad92 363void AliCFGrid::SetElementError(Double_t *var, Float_t val)
563113d0 364{
365 //
366 // Set squared error to val in a bin corresponding to a set of input variables
367 //
368 Int_t unfl=0;
369 Int_t ovfl=0;
370 for(Int_t i=0;i<fNVar;i++){
371 Int_t nbins=fNVarBins[i]+1;
372 Float_t *bins=new Float_t[nbins];
373 for(Int_t ibin =0;ibin<nbins;ibin++){
374 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
375 }
376
377 fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
378 //underflows
379
380 if(var[i] < bins[0]){
381 unfl=1;
382 }
383
384 //overflows
385
386 if(var[i] > bins[nbins-1]){
387 ovfl=1;
388 }
389 delete [] bins;
390 }
391
1e9dad92 392 //move to the TH/THnSparse convention in N-dim bin numbering
393 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
394
563113d0 395 if(!(ovfl==1 || unfl==1)){
396 SetElementError(fIndex,val);
397 }
398 else{
399 AliInfo(Form(" input variables outside the grid, no value set"));
400 }
401}
402//____________________________________________________________________
1e9dad92 403void AliCFGrid::Fill(Double_t *var, Double_t weight)
563113d0 404{
1e9dad92 405
563113d0 406 //
1e9dad92 407 // Fill the grid,
408 // given a set of values of the input variable,
409 // with weight (by default w=1)
563113d0 410 //
1e9dad92 411
412
413 Int_t isunfl=0;
414 Int_t isovfl=0;
415 Int_t *unfl=new Int_t[fNVar];
416 Int_t *ovfl=new Int_t[fNVar];
417
418 for(Int_t i=0;i<fNVar;i++){
419 unfl[i]=0;
420 ovfl[i]=0;
563113d0 421 }
1e9dad92 422
563113d0 423 for(Int_t i=0;i<fNVar;i++){
424 Int_t nbins=fNVarBins[i]+1;
425 Float_t *bins=new Float_t[nbins];
426 for(Int_t ibin =0;ibin<nbins;ibin++){
427 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
428 }
429
430 fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
431 //underflows
432
433 if(var[i] < bins[0]){
1e9dad92 434 unfl[i]=1;
435 isunfl=1;
563113d0 436 }
437
438 //overflows
439
440 if(var[i] > bins[nbins-1]){
1e9dad92 441 ovfl[i]=1;
442 isovfl=1;
563113d0 443 }
444 delete [] bins;
445 }
446
1e9dad92 447 //exclusive under/overflows
563113d0 448
563113d0 449 for(Int_t i=0;i<fNVar;i++){
1e9dad92 450 Bool_t add=kTRUE;
451 for(Int_t j=0;j<fNVar;j++){
452 if(i==j)continue;
453 if(!(unfl[j]==0 && ovfl[j]==0))add=kFALSE;
563113d0 454 }
1e9dad92 455 if(add && unfl[i]==1)fNunfl[i]++;
456 if(add && ovfl[i]==1)fNovfl[i]++;
563113d0 457 }
458
1e9dad92 459 delete [] unfl;
460 delete [] ovfl;
461
563113d0 462 // Total number of entries, overflows and underflows
463
464 fNentriesTot++;
1e9dad92 465 if(isunfl)fNunflTot++;
466 if(isovfl)fNovflTot++;
563113d0 467
468 //if not ovfl/unfl, fill the element
1e9dad92 469
470 if(!(isovfl==1 || isunfl==1)){
563113d0 471 Int_t ind =GetBinIndex(fIndex);
472 fData[ind]+=weight;
473 if(fSumW2)fErr2[ind]+=(weight*weight);
474 }
475}
476//____________________________________________________________________
477Float_t AliCFGrid::GetOverFlows( Int_t ivar) const {
478 //
479 // Get overflows in variable var
480 //
481 return fNovfl[ivar];
482}
483//____________________________________________________________________
484Float_t AliCFGrid::GetOverFlows() const {
485 //
486 // Get overflows
487 //
488 return fNovflTot;
489}
490//____________________________________________________________________
491Float_t AliCFGrid::GetUnderFlows( Int_t ivar) const {
492 //
493 // Get overflows in variable var
494 //
495 return fNunfl[ivar];
496}
497//____________________________________________________________________
498Float_t AliCFGrid::GetUnderFlows() const {
499 //
500 // Get overflows
501 //
502 return fNunflTot;
503}
504//____________________________________________________________________
505Float_t AliCFGrid::GetEntries( ) const {
506 //
507 // Get total entries (in grid + over/underflows)
508 //
509 return fNentriesTot;
510}
563113d0 511//___________________________________________________________________
1e9dad92 512TH1D *AliCFGrid::Project(Int_t ivar) const
563113d0 513{
514 //
515 // Make a 1D projection along variable ivar
516
517
518 Int_t nbins =fNVarBins[ivar];
519 Float_t *bins = new Float_t[nbins+1];
520 for (Int_t i=0;i<=fNVar;i++){
521 }
522 for(Int_t ibin =0;ibin<nbins+1;ibin++){
523 bins[ibin] = fVarBinLimits[ibin+fOffset[ivar]];
524 }
525
526 char pname[40];
527 sprintf(pname,"%s%s_%i",GetName(),"_proj1D_var", ivar);
528 char htitle[40];
529 sprintf(htitle,"%s%s_%i",GetName(),"_proj1D_var", ivar);
530
1e9dad92 531 TH1D *proj1D=0;
563113d0 532
533 //check if a projection with identical name exist
534 TObject *obj = gROOT->FindObject(pname);
1e9dad92 535 if (obj && obj->InheritsFrom("TH1D")) {
536 proj1D = (TH1D*)obj;
563113d0 537 proj1D->Reset();
538 }
539
540 if(!proj1D){
1e9dad92 541 proj1D =new TH1D(pname,htitle, nbins, bins);
563113d0 542 }
543
544 delete [] bins;
545 Float_t sum=0;
546 Float_t *data= new Float_t[nbins];
547 Float_t *err= new Float_t[nbins];
548
549 for(Int_t ibin=0;ibin<nbins;ibin++)data[ibin]=0;
550 for(Int_t ibin=0;ibin<nbins;ibin++)err[ibin]=0;
551 for(Int_t iel=0;iel<fNDim;iel++){
552 data[GetBinIndex(ivar,iel)]+=fData[iel];
553 if(fSumW2)err[GetBinIndex(ivar,iel)]+=fErr2[iel];
554 }
555
556 for(Int_t ibin =0;ibin<nbins;ibin++){
557 proj1D->SetBinContent(ibin+1,data[ibin]);
558 proj1D->SetBinError(ibin+1,TMath::Sqrt(data[ibin]));
559 if(fSumW2)proj1D->SetBinError(ibin+1,TMath::Sqrt(err[ibin]));
560 sum+=data[ibin];
561 }
562
563 delete [] data;
564 delete [] err;
565 proj1D->SetBinContent(nbins+1,GetOverFlows(ivar));
566 proj1D->SetBinContent(0,GetUnderFlows(ivar));
1e9dad92 567 proj1D->SetEntries(sum+GetUnderFlows(ivar)+GetOverFlows(ivar));
563113d0 568 return proj1D;
569}
570
571//___________________________________________________________________
1e9dad92 572TH2D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2) const
563113d0 573{
574 //
575 // Make a 2D projection along variable ivar
576
577 Int_t nbins1 =fNVarBins[ivar1];
578 Int_t nbins2 =fNVarBins[ivar2];
579
580 Float_t *bins1 = new Float_t[nbins1+1];
581 Float_t *bins2 = new Float_t[nbins2+1];
582
583 for(Int_t ibin =0;ibin<nbins1+1;ibin++){
584 bins1[ibin] = fVarBinLimits[ibin+fOffset[ivar1]];
585 }
586 for(Int_t ibin =0;ibin<nbins2+1;ibin++){
587 bins2[ibin] = fVarBinLimits[ibin+fOffset[ivar2]];
588 }
589
590 char pname[40];
591 sprintf(pname,"%s%s_%i_%i",GetName(),"_proj2D_var",ivar1,ivar2);
592 char htitle[40];
593 sprintf(htitle,"%s%s_%i_%i",GetName(),"_proj2D_var",ivar1,ivar2);
594
1e9dad92 595 TH2D *proj2D=0;
563113d0 596
597 //check if a projection with identical name exist
598 TObject *obj = gROOT->FindObject(pname);
1e9dad92 599 if (obj && obj->InheritsFrom("TH2D")) {
600 proj2D = (TH2D*)obj;
563113d0 601 proj2D->Reset();
602 }
603
604 if(!proj2D){
1e9dad92 605 proj2D =new TH2D(pname,htitle, nbins1, bins1,nbins2,bins2);
563113d0 606 }
607
608 delete [] bins1;
609 delete [] bins2;
610
611
612 Float_t sum=0;
613 Float_t **data=new Float_t*[nbins1];
614 Float_t *data2=new Float_t[nbins1*nbins2];
615 Float_t **err=new Float_t*[nbins1];
616 Float_t *err2=new Float_t[nbins1*nbins2];
617 for(Int_t i=0;i<nbins1;i++)data[i] = data2+i*nbins2;
618 for(Int_t i=0;i<nbins1;i++)err[i] = err2+i*nbins2;
619
620 for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
621 for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
622 data[ibin1][ibin2]=0;
623 err[ibin1][ibin2]=0;
624 }
625 }
626
627 for(Int_t iel=0;iel<fNDim;iel++){
628 data[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)]+=fData[iel];
629 if(fSumW2)err[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)]+=fErr2[iel];
630 }
631
632 for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
633 for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
634 proj2D->SetBinContent(ibin1+1,ibin2+1,data[ibin1][ibin2]);
635 proj2D->SetBinError(ibin1+1,ibin2+1,TMath::Sqrt(data[ibin1][ibin2]));
636 if(fSumW2)proj2D->SetBinError(ibin1+1,ibin2+1,TMath::Sqrt(err[ibin1][ibin2]));
637 sum+=data[ibin1][ibin2];
638 }
639
640 }
641 delete data;
642 delete data2;
643 delete err;
644 delete err2;
645 proj2D->SetBinContent(0,nbins2/2,GetUnderFlows(ivar1));
646 proj2D->SetBinContent(nbins1+1,nbins2/2,GetOverFlows(ivar1));
647 proj2D->SetBinContent(nbins1/2,0,GetUnderFlows(ivar2));
648 proj2D->SetBinContent(nbins1/2,nbins2+1,GetOverFlows(ivar2));
1e9dad92 649 proj2D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2));
563113d0 650 return proj2D;
651}
652//___________________________________________________________________
1e9dad92 653TH3D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
563113d0 654{
655 //
656 // Make a 3D projection along variable ivar
657
658 Int_t nbins1 =fNVarBins[ivar1];
659 Int_t nbins2 =fNVarBins[ivar2];
660 Int_t nbins3 =fNVarBins[ivar3];
661
662 Float_t *bins1 = new Float_t[nbins1+1];
663 Float_t *bins2 = new Float_t[nbins2+1];
664 Float_t *bins3 = new Float_t[nbins3+1];
665
666 for(Int_t ibin =0;ibin<nbins1+1;ibin++){
667 bins1[ibin] = fVarBinLimits[ibin+fOffset[ivar1]];
668 }
669 for(Int_t ibin =0;ibin<nbins2+1;ibin++){
670 bins2[ibin] = fVarBinLimits[ibin+fOffset[ivar2]];
671 }
672 for(Int_t ibin =0;ibin<nbins3+1;ibin++){
673 bins3[ibin] = fVarBinLimits[ibin+fOffset[ivar3]];
674 }
675
676 char pname[40];
677 sprintf(pname,"%s%s_%i_%i_%i",GetName(),"_proj3D_var",ivar1,ivar2,ivar3);
678 char htitle[40];
679 sprintf(htitle,"%s%s_%i_%i_%i",GetName(),"_proj3D_var",ivar1,ivar2,ivar3);
680
1e9dad92 681 TH3D *proj3D=0;
563113d0 682
683 //check if a projection with identical name exist
684 TObject *obj = gROOT->FindObject(pname);
1e9dad92 685 if (obj && obj->InheritsFrom("TH3D")) {
686 proj3D = (TH3D*)obj;
563113d0 687 proj3D->Reset();
688 }
689
690 if(!proj3D){
1e9dad92 691 proj3D =new TH3D(pname,htitle, nbins1,bins1,nbins2,bins2,nbins3,bins3);
563113d0 692 }
693
694 delete [] bins1;
695 delete [] bins2;
696 delete [] bins3;
697
698
699 Float_t sum=0;
700 Float_t ***data=new Float_t**[nbins1];
701 Float_t **data2=new Float_t*[nbins1*nbins2];
702 Float_t *data3=new Float_t[nbins1*nbins2*nbins3];
703 Float_t ***err=new Float_t**[nbins1];
704 Float_t **err2=new Float_t*[nbins1*nbins2];
705 Float_t *err3=new Float_t[nbins1*nbins2*nbins3];
706 for(Int_t i=0;i<nbins1;i++)data[i] = data2+i*nbins2;
707 for(Int_t i=0;i<nbins1;i++)err[i] = err2+i*nbins2;
708 for(Int_t i=0;i<nbins1;i++){
709 for(Int_t j=0;j<nbins2;j++){
710 data[i][j] = data3+i*nbins2*nbins3+j*nbins3;
711 err[i][j] = err3+i*nbins2*nbins3+j*nbins3;
712 }
713 }
714 for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
715 for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
716 for(Int_t ibin3 =0;ibin3<nbins3;ibin3++){
717 data[ibin1][ibin2][ibin3]=0;
718 err[ibin1][ibin2][ibin3]=0;
719 }
720 }
721 }
722
723 for(Int_t iel=0;iel<fNDim;iel++){
724 data[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)][GetBinIndex(ivar3,iel)]+=fData[iel];
725 if(fSumW2)err[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)][GetBinIndex(ivar3,iel)]+=fErr2[iel];
726 }
727
728 for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
729 for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
730 for(Int_t ibin3 =0;ibin3<nbins3;ibin3++){
731 proj3D->SetBinContent(ibin1+1,ibin2+1,ibin3+1,data[ibin1][ibin2][ibin3]);
732 proj3D->SetBinError(ibin1+1,ibin2+1,ibin3+1,TMath::Sqrt(data[ibin1][ibin2][ibin3]));
733 if(fSumW2)proj3D->SetBinError(ibin1+1,ibin2+1,ibin3+1,TMath::Sqrt(err[ibin1][ibin2][ibin3]));
734 sum+=data[ibin1][ibin2][ibin3];
735 }
736 }
737 }
738
739 delete data;
740 delete data2;
741 delete data3;
742 delete err;
743 delete err2;
744 delete err3;
1e9dad92 745 proj3D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2)+GetUnderFlows(ivar3)+GetOverFlows(ivar3));
563113d0 746 return proj3D;
747}
748
749//___________________________________________________________________
1e9dad92 750TH1D *AliCFGrid::Slice(Int_t ivar, Double_t *varMin, Double_t* varMax) const
563113d0 751{
752 //
753 // Make a slice along variable ivar in range [varMin,varMax]
754
755
756 Int_t nbins =fNVarBins[ivar];
757 Float_t *bins = new Float_t[nbins+1];
758 for (Int_t i=0;i<=fNVar;i++){
759 }
760 for(Int_t ibin =0;ibin<nbins+1;ibin++){
761 bins[ibin] = fVarBinLimits[ibin+fOffset[ivar]];
762 }
763
764 char pname[40];
765 sprintf(pname,"%s%s_%i",GetName(),"_proj1D_var", ivar);
766 char htitle[40];
767 sprintf(htitle,"%s%s_%i",GetName(),"_proj1D_var", ivar);
768
1e9dad92 769 TH1D *proj1D=0;
563113d0 770
771 //check if a projection with identical name exist
772 TObject *obj = gROOT->FindObject(pname);
1e9dad92 773 if (obj && obj->InheritsFrom("TH1D")) {
774 proj1D = (TH1D*)obj;
563113d0 775 proj1D->Reset();
776 }
777
778 if(!proj1D){
1e9dad92 779 proj1D =new TH1D(pname,htitle, nbins, bins);
563113d0 780 }
781
782 delete [] bins;
783
784
785 Int_t *indexMin=new Int_t[fNVar];
786 Int_t *indexMax=new Int_t[fNVar];
787
788
789 //Find out the min and max bins
790
791 for(Int_t i=0;i<fNVar;i++){
792 Float_t xmin=varMin[i]; // the min values
793 Float_t xmax=varMax[i]; // the min values
794 Int_t nbins=fNVarBins[i]+1;
795 Float_t *bins=new Float_t[nbins];
796 for(Int_t ibin =0;ibin<nbins;ibin++){
797 bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
798 }
799 indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
800 indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
563113d0 801 delete [] bins;
802 }
803
804 Float_t sum=0;
805 Float_t *data= new Float_t[nbins];
806 for(Int_t ibin=0;ibin<nbins;ibin++)data[ibin]=0;
807
808 Int_t *index= new Int_t[fNVar];
809 Int_t ielmin=GetBinIndex(indexMin);
810 Int_t ielmax=GetBinIndex(indexMax);
811 for(Int_t iel=ielmin;iel<=ielmax;iel++){
812 GetBinIndex(iel,index);
813 Bool_t isIn=kTRUE;
814 for (Int_t j=0;j<fNVar;j++){
815 if(!(index[j]>=indexMin[j] && index[j]<=indexMax[j]))isIn=kFALSE;
816 break;
817 }
818 if(isIn)data[GetBinIndex(ivar,iel)]+=fData[iel];
819 }
820
821 delete [] index;
822
823
824 for(Int_t ibin =0;ibin<nbins;ibin++){
825 proj1D->SetBinContent(ibin+1,data[ibin]);
826 proj1D->SetBinError(ibin+1,TMath::Sqrt(data[ibin]));
827 sum+=data[ibin];
828 }
829
830 delete [] data;
831
832 proj1D->SetEntries(sum);
833 return proj1D;
834}
835
563113d0 836
837//____________________________________________________________________
1e9dad92 838void AliCFGrid::Add(AliCFVGrid* aGrid, Double_t c)
563113d0 839{
840 //
841 //add aGrid to the current one
842 //
843
844 if(aGrid->GetNVar()!=fNVar){
845 AliInfo("Different number of variables, cannot add the grids");
846 return;
847 }
848 if(aGrid->GetNDim()!=fNDim){
849 AliInfo("Different number of dimensions, cannot add the grids!");
850 return;
851 }
852
853 if(!fSumW2 && aGrid->GetSumW2())SumW2();
854
855 for(Int_t iel=0;iel<fNDim;iel++){
856 fData[iel]+=(c*aGrid->GetElement(iel));
857 if(fSumW2){
1e9dad92 858 Float_t err=aGrid->GetElementError(iel);
563113d0 859 fErr2[iel]+=c*c*err*err;
860 }
861 }
862
863 //Add entries, overflows and underflows
864
865 fNentriesTot+= c*aGrid->GetEntries();
866 fNunflTot+= c*aGrid->GetUnderFlows();
867 fNovflTot+= c*aGrid->GetOverFlows();
868 for(Int_t j=0;j<fNVar;j++){
869 fNunfl[j]+= c*aGrid->GetUnderFlows(j);
870 fNovfl[j]+= c*aGrid->GetUnderFlows(j);
871 }
872}
873//____________________________________________________________________
1e9dad92 874void AliCFGrid::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
563113d0 875{
876 //
877 //add aGrid1 and aGrid2
878 //
879
880 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
881 AliInfo("Different number of variables, cannot add the grids");
882 return;
883 }
884 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
885 AliInfo("Different number of dimensions, cannot add the grids!");
886 return;
887 }
888
889 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
890
891 Float_t cont1,cont2,err1,err2;
892
893 for(Int_t iel=0;iel<fNDim;iel++){
894 cont1=aGrid1->GetElement(iel);
895 cont2=aGrid2->GetElement(iel);
896 SetElement(iel,c1*cont1+c2*cont2);
897 if(fSumW2){
1e9dad92 898 err1=aGrid1->GetElementError(iel);
899 err2=aGrid2->GetElementError(iel);
900 SetElementError(iel,TMath::Sqrt(c1*c1*err1*err1+c2*c2*err2*err2));
563113d0 901 }
902 }
903
904 //Add entries, overflows and underflows
905
906 fNentriesTot= c1*aGrid1->GetEntries()+c2*aGrid2->GetEntries();
907 fNunflTot= c1*aGrid1->GetUnderFlows()+c2*aGrid2->GetUnderFlows();
908 fNovflTot= c1*aGrid1->GetOverFlows()+c2*aGrid2->GetOverFlows();
909 for(Int_t j=0;j<fNVar;j++){
910 fNunfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
911 fNovfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
912 }
913}
914//____________________________________________________________________
1e9dad92 915void AliCFGrid::Multiply(AliCFVGrid* aGrid, Double_t c)
563113d0 916{
917 //
918 //multiply grid aGrid by the current one
919 //
920
921 if(aGrid->GetNVar()!=fNVar){
922 AliInfo("Different number of variables, cannot multiply the grids");
923 return;
924 }
925 if(aGrid->GetNDim()!=fNDim){
926 AliInfo("Different number of dimensions, cannot multiply the grids!");
927 return;
928 }
929
930 if(!fSumW2 && aGrid->GetSumW2())SumW2();
931
932 Float_t cont1,cont2,err1,err2;
933
934 for(Int_t iel=0;iel<fNDim;iel++){
935 cont1=GetElement(iel);
936 cont2=c*aGrid->GetElement(iel);
937 SetElement(iel,cont1*cont2);
938 if(fSumW2){
1e9dad92 939 err1=GetElementError(iel);
940 err2=aGrid->GetElementError(iel);
941 SetElementError(iel,TMath::Sqrt(c*c*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
563113d0 942 }
943 }
944
945 //Set entries to the number of bins, preserve original overflows and underflows
946
947 fNentriesTot=fNDim;
948 fNunflTot=GetUnderFlows();
949 fNovflTot=GetOverFlows();
950 for(Int_t j=0;j<fNVar;j++){
951 fNunfl[j]= GetUnderFlows(j);
952 fNovfl[j]= GetUnderFlows(j);
953 }
954}
955//____________________________________________________________________
1e9dad92 956void AliCFGrid::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1, Double_t c2)
563113d0 957{
958 //
959 //multiply grids aGrid1 and aGrid2
960 //
961
962 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
963 AliInfo("Different number of variables, cannot multiply the grids");
964 return;
965 }
966 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
967 AliInfo("Different number of dimensions, cannot multiply the grids!");
968 return;
969 }
970
971 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
972
973 Float_t cont1,cont2,err1,err2;
974 for(Int_t iel=0;iel<fNDim;iel++){
975 cont1=c1*aGrid1->GetElement(iel);
976 cont2=c2*aGrid2->GetElement(iel);
977 SetElement(iel,cont1*cont2);
978 if(fSumW2){
1e9dad92 979 err1=aGrid1->GetElementError(iel);
980 err2=aGrid2->GetElementError(iel);
981 SetElementError(iel,TMath::Sqrt(c1*c1*c2*c2*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
563113d0 982 }
983 }
984
985 //Set entries to the number of bins, preserve original overflows and underflows
986
987 fNentriesTot=fNDim;
988 fNunflTot=GetUnderFlows();
989 fNovflTot=GetOverFlows();
990 for(Int_t j=0;j<fNVar;j++){
991 fNunfl[j]= GetUnderFlows(j);
992 fNovfl[j]= GetUnderFlows(j);
993 }
994}
995//____________________________________________________________________
1e9dad92 996void AliCFGrid::Divide(AliCFVGrid* aGrid, Double_t c)
563113d0 997{
998 //
999 //divide current grid by grid aGrid
1000 //
1001
563113d0 1002 if(aGrid->GetNVar()!=fNVar){
1003 AliInfo("Different number of variables, cannot divide the grids");
1004 return;
1005 }
1006 if(aGrid->GetNDim()!=fNDim){
1007 AliInfo("Different number of dimensions, cannot divide the grids!");
1008 return;
1009 }
1010 if(!c){AliInfo(Form("c is %f, cannot divide!",c)); return;}
1011
1012 if(!fSumW2 && aGrid->GetSumW2())SumW2();
1013
1e9dad92 1014 Float_t cont1,cont2,err1,err2,den;
563113d0 1015 for(Int_t iel=0;iel<fNDim;iel++){
1016 cont1=GetElement(iel);
1017 cont2=aGrid->GetElement(iel);
1018 if(cont2)SetElement(iel,cont1/(c*cont2));
1019 else SetElement(iel,0);
1020 if(fSumW2){
1e9dad92 1021 err1=GetElementError(iel);
1022 err2=aGrid->GetElementError(iel);
563113d0 1023 if(!cont2){SetElementError(iel,0.); continue;}
1e9dad92 1024 den=cont2*cont2*cont2*c*c;
1025 SetElementError(iel,TMath::Sqrt((cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
563113d0 1026 }
1027 }
1e9dad92 1028
563113d0 1029 //Set entries to the number of bins, preserve original overflows and underflows
1030
1031 fNentriesTot=fNDim;
1032 fNunflTot=GetUnderFlows();
1033 fNovflTot=GetOverFlows();
1034 for(Int_t j=0;j<fNVar;j++){
1035 fNunfl[j]= GetUnderFlows(j);
1036 fNovfl[j]= GetUnderFlows(j);
1037 }
1038}
1039//____________________________________________________________________
1e9dad92 1040void AliCFGrid::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
563113d0 1041{
1042 //
1043 //divide grids aGrid1,aGrid2
1044 //
1045
1046 TString opt = option;
1047 opt.ToUpper();
1048
1049 if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
1050 AliInfo("Different number of variables, cannot divide the grids");
1051 return;
1052 }
1053 if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
1054 AliInfo("Different number of dimensions, cannot divide the grids!");
1055 return;
1056 }
1057 if(!c2){AliInfo(Form("c2 is %f, cannot divide!",c2)); return;}
1058
1059 if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
1060
1061 Float_t cont1,cont2,err1,err2,r,den;
1062
1063 for(Int_t iel=0;iel<fNDim;iel++){
1064 cont1=aGrid1->GetElement(iel);
1065 cont2=aGrid2->GetElement(iel);
1066 if(cont2)SetElement(iel,c1*cont1/(c2*cont2));
1067 else SetElement(iel,0);
1068 if(fSumW2){
1e9dad92 1069 err1=aGrid1->GetElementError(iel);
1070 err2=aGrid2->GetElementError(iel);
563113d0 1071 if(!cont2){SetElementError(iel,0.); continue;}
1072 if (opt.Contains("B")){
1073 if(cont1!=cont2){
1074 r=cont1/cont2;
1e9dad92 1075 SetElementError(iel,TMath::Sqrt(TMath::Abs(((1.-2.*r)*err1*err1+r*r*err2*err2)/(cont2*cont2))));
563113d0 1076 }else{
1077 SetElementError(iel,0.);
1078 }
1079 }else{
1080 den=cont2*cont2*cont2*cont2*c2*c2;
1e9dad92 1081 SetElementError(iel,TMath::Sqrt(c1*c1*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
563113d0 1082 }
1083 }
1084 }
1085
1086 //Set entries to the number of bins, preserve original overflows and underflows
1087
1088 fNentriesTot=fNDim;
1089 fNunflTot=GetUnderFlows();
1090 fNovflTot=GetOverFlows();
1091 for(Int_t j=0;j<fNVar;j++){
1092 fNunfl[j]= GetUnderFlows(j);
1093 fNovfl[j]= GetUnderFlows(j);
1094 }
1095}
1096//____________________________________________________________________
1097void AliCFGrid::SumW2()
1098{
1099 //
1100 //set calculation of the squared sum of the weighted entries
1101 //
1102 if(!fSumW2){
1103 fErr2=new Float_t [fNDim];
1104 //init....
1105 for(Int_t iel=0;iel<fNDim;iel++){
1106 fErr2[iel]=fData[iel];
1107 }
1108 }
1109
1110 fSumW2=kTRUE;
1111}
563113d0 1112//____________________________________________________________________
1113void AliCFGrid::Copy(TObject& c) const
1114{
1115 //
1116 // copy function
1117 //
1118 AliCFGrid& target = (AliCFGrid &) c;
1119
563113d0 1120 target.fNunflTot = fNunflTot;
1121 target.fNovflTot = fNovflTot;
1122 target.fNentriesTot = fNentriesTot;
563113d0 1123 if (fNunfl)
1124 target.fNunfl = fNunfl;
1125 if (fNovfl)
1126 target.fNovfl = fNovfl;
563113d0 1127 if (fData)
1128 target.fData = fData;
1129 if (fErr2)
1130 target.fErr2 = fErr2;
1131
1132}