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