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