New support for QA histos
[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 eventually replaced by AliCFGridSparse  //
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     Double_t *bins=new Double_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     Double_t *bins=new Double_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     Double_t *bins=new Double_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     Double_t *bins=new Double_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     Double_t *bins=new Double_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     Double_t xmin=varMin[i]; // the min values  
771     Double_t xmax=varMax[i]; // the max values  
772     Int_t nbins=fNVarBins[i]+1;
773     Double_t *bins=new Double_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
790
791   for(Int_t iel=ielmin;iel<=ielmax;iel++){
792     GetBinIndex(iel,index);
793     Bool_t isIn=kTRUE;
794     for (Int_t j=0;j<fNVar;j++){
795       if(!(index[j]>=indexMin[j] && index[j]<=indexMax[j]))isIn=kFALSE;   
796       break;
797     }
798     if(isIn)data[GetBinIndex(ivar,iel)]+=fData[iel];
799   }
800
801   delete [] index;
802   delete [] indexMin;
803   delete [] indexMax;
804
805
806   for(Int_t ibin =0;ibin<nbins;ibin++){
807     proj1D->SetBinContent(ibin+1,data[ibin]);
808     proj1D->SetBinError(ibin+1,TMath::Sqrt(data[ibin]));
809     sum+=data[ibin];
810   }
811
812   delete [] data;
813
814   proj1D->SetEntries(sum);
815   return proj1D;
816
817
818
819 //____________________________________________________________________
820 void AliCFGrid::Add(AliCFVGrid* aGrid, Double_t c)
821 {
822   //
823   //add aGrid to the current one
824   //
825
826   if(aGrid->GetNVar()!=fNVar){
827     AliInfo("Different number of variables, cannot add the grids");
828     return;
829   } 
830   if(aGrid->GetNDim()!=fNDim){
831     AliInfo("Different number of dimensions, cannot add the grids!");
832     return;
833   } 
834   
835   if(!fSumW2  && aGrid->GetSumW2())SumW2();
836
837   for(Int_t iel=0;iel<fNDim;iel++){
838     fData[iel]+=(c*aGrid->GetElement(iel));
839     if(fSumW2){
840       Float_t err=aGrid->GetElementError(iel);  
841       fErr2[iel]+=c*c*err*err;
842     }
843   }
844
845   //Add entries, overflows and underflows
846
847   fNentriesTot+= c*aGrid->GetEntries();
848   for(Int_t j=0;j<fNVar;j++){
849     fNunfl[j]+= c*aGrid->GetUnderFlows(j);
850     fNovfl[j]+= c*aGrid->GetOverFlows(j);
851   }
852 }
853 //____________________________________________________________________
854 void AliCFGrid::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
855 {
856   //
857   //add aGrid1 and aGrid2
858   //
859
860   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
861     AliInfo("Different number of variables, cannot add the grids");
862     return;
863   } 
864   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
865     AliInfo("Different number of dimensions, cannot add the grids!");
866     return;
867   } 
868   
869   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
870
871   Float_t cont1,cont2,err1,err2;  
872
873   for(Int_t iel=0;iel<fNDim;iel++){
874     cont1=aGrid1->GetElement(iel);
875     cont2=aGrid2->GetElement(iel);
876     SetElement(iel,c1*cont1+c2*cont2);
877     if(fSumW2){
878       err1=aGrid1->GetElementError(iel);
879       err2=aGrid2->GetElementError(iel);
880       SetElementError(iel,TMath::Sqrt(c1*c1*err1*err1+c2*c2*err2*err2));
881     }
882   }
883
884   //Add entries, overflows and underflows
885
886   fNentriesTot= c1*aGrid1->GetEntries()+c2*aGrid2->GetEntries();
887   for(Int_t j=0;j<fNVar;j++){
888     fNunfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
889     fNovfl[j]= c1*aGrid1->GetOverFlows(j)+c2*aGrid2->GetOverFlows(j);
890   }
891 }
892 //____________________________________________________________________
893 void AliCFGrid::Multiply(AliCFVGrid* aGrid, Double_t c)
894 {
895   //
896   //multiply grid aGrid by the current one
897   //
898
899   if(aGrid->GetNVar()!=fNVar){
900     AliInfo("Different number of variables, cannot multiply the grids");
901     return;
902   } 
903   if(aGrid->GetNDim()!=fNDim){
904     AliInfo("Different number of dimensions, cannot multiply the grids!");
905     return;
906   } 
907   
908   if(!fSumW2  && aGrid->GetSumW2())SumW2();
909   
910   Float_t cont1,cont2,err1,err2;  
911
912   for(Int_t iel=0;iel<fNDim;iel++){
913     cont1=GetElement(iel);  
914     cont2=c*aGrid->GetElement(iel);  
915     SetElement(iel,cont1*cont2);
916     if(fSumW2){
917       err1=GetElementError(iel);  
918       err2=aGrid->GetElementError(iel);  
919       SetElementError(iel,TMath::Sqrt(c*c*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
920     }
921   }
922
923   //Set entries to the number of bins, preserve original overflows and underflows
924
925   fNentriesTot=fNDim;
926   for(Int_t j=0;j<fNVar;j++){
927     fNunfl[j]= GetUnderFlows(j);
928     fNovfl[j]= GetOverFlows(j);
929   }
930 }
931 //____________________________________________________________________
932 void AliCFGrid::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1, Double_t c2)
933 {
934   //
935   //multiply grids aGrid1 and aGrid2
936   //
937
938   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
939     AliInfo("Different number of variables, cannot multiply the grids");
940     return;
941   } 
942   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
943     AliInfo("Different number of dimensions, cannot multiply the grids!");
944     return;
945   } 
946   
947   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
948
949   Float_t cont1,cont2,err1,err2;  
950   for(Int_t iel=0;iel<fNDim;iel++){
951     cont1=c1*aGrid1->GetElement(iel);  
952     cont2=c2*aGrid2->GetElement(iel);  
953     SetElement(iel,cont1*cont2);
954     if(fSumW2){
955       err1=aGrid1->GetElementError(iel);  
956       err2=aGrid2->GetElementError(iel);  
957       SetElementError(iel,TMath::Sqrt(c1*c1*c2*c2*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
958     }
959   }
960
961   //Set entries to the number of bins, preserve original overflows and underflows
962
963   fNentriesTot=fNDim;
964   for(Int_t j=0;j<fNVar;j++){
965     fNunfl[j]= GetUnderFlows(j);
966     fNovfl[j]= GetOverFlows(j);
967   }
968 }
969 //____________________________________________________________________
970 void AliCFGrid::Divide(AliCFVGrid* aGrid, Double_t c)
971 {
972   //
973   //divide current grid by grid aGrid
974   //
975
976   if(aGrid->GetNVar()!=fNVar){
977     AliInfo("Different number of variables, cannot divide the grids");
978     return;
979   } 
980   if(aGrid->GetNDim()!=fNDim){
981     AliInfo("Different number of dimensions, cannot divide the grids!");
982     return;
983   } 
984  if(!c){AliInfo(Form("c is %f, cannot divide!",c)); return;} 
985   
986   if(!fSumW2  && aGrid->GetSumW2())SumW2();
987
988   Float_t cont1,cont2,err1,err2,den;
989   for(Int_t iel=0;iel<fNDim;iel++){
990     cont1=GetElement(iel);  
991     cont2=aGrid->GetElement(iel);
992     if(cont2)SetElement(iel,cont1/(c*cont2));
993     else SetElement(iel,0);
994     if(fSumW2){
995       err1=GetElementError(iel);  
996       err2=aGrid->GetElementError(iel);  
997       if(!cont2){SetElementError(iel,0.); continue;}
998       den=cont2*cont2*cont2*c*c;
999       SetElementError(iel,TMath::Sqrt((cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
1000     }
1001   }
1002   
1003   //Set entries to the number of bins, preserve original overflows and underflows
1004
1005   fNentriesTot=fNDim;
1006   for(Int_t j=0;j<fNVar;j++){
1007     fNunfl[j]= GetUnderFlows(j);
1008     fNovfl[j]= GetOverFlows(j);
1009   }
1010 }
1011 //____________________________________________________________________
1012 void AliCFGrid::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
1013 {
1014   //
1015   //divide grids aGrid1,aGrid2
1016   //
1017
1018   TString opt = option;
1019   opt.ToUpper();
1020
1021   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
1022     AliInfo("Different number of variables, cannot divide the grids");
1023     return;
1024   }
1025   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
1026     AliInfo("Different number of dimensions, cannot divide the grids!");
1027     return;
1028   }
1029   if(!c2){AliInfo(Form("c2 is %f, cannot divide!",c2)); return;} 
1030   
1031   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
1032
1033   Float_t cont1,cont2,err1,err2,r,den;
1034  
1035   for(Int_t iel=0;iel<fNDim;iel++){
1036     cont1=aGrid1->GetElement(iel);  
1037     cont2=aGrid2->GetElement(iel);  
1038     if(cont2)SetElement(iel,c1*cont1/(c2*cont2));
1039     else SetElement(iel,0);
1040      if(fSumW2){
1041       err1=aGrid1->GetElementError(iel);  
1042       err2=aGrid2->GetElementError(iel);  
1043       if(!cont2){SetElementError(iel,0.); continue;}
1044       if (opt.Contains("B")){
1045         if(cont1!=cont2){
1046           r=cont1/cont2;            
1047           SetElementError(iel,TMath::Sqrt(TMath::Abs(((1.-2.*r)*err1*err1+r*r*err2*err2)/(cont2*cont2))));
1048         }else{
1049           SetElementError(iel,0.);
1050         }
1051       }else{
1052         den=cont2*cont2*cont2*cont2*c2*c2;
1053         SetElementError(iel,TMath::Sqrt(c1*c1*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
1054       }
1055     }
1056   }
1057
1058   //Set entries to the number of bins, preserve original overflows and underflows
1059
1060   fNentriesTot=fNDim;
1061   for(Int_t j=0;j<fNVar;j++){
1062     fNunfl[j]= GetUnderFlows(j);
1063     fNovfl[j]= GetOverFlows(j);
1064   }
1065 }
1066 //____________________________________________________________________
1067 void AliCFGrid::Rebin(const Int_t* group)
1068 {
1069   //
1070   // Not yet implemented
1071   //
1072   for(Int_t i=0;i<fNVar;i++){
1073     if(group[i]!=1)AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
1074   }
1075   AliInfo(Form("This method was so far not implemented for AliCFGrid, but it is available for AliCFGridSparse"));
1076
1077 }
1078 //____________________________________________________________________
1079 void AliCFGrid::SumW2()
1080 {
1081   //
1082   //set calculation of the squared sum of the weighted entries
1083   //
1084   if(!fSumW2){
1085     fErr2=new Float_t [fNDim];
1086     //init....
1087     for(Int_t iel=0;iel<fNDim;iel++){
1088       fErr2[iel]=fData[iel];
1089     }
1090   }
1091
1092   fSumW2=kTRUE;
1093 }
1094 //____________________________________________________________________
1095 void AliCFGrid::Copy(TObject& c) const
1096 {
1097   //
1098   // copy function
1099   //
1100   AliCFGrid& target = (AliCFGrid &) c;
1101
1102   target.fNentriesTot = fNentriesTot;
1103   if (fNunfl)
1104     target.fNunfl = fNunfl;
1105   if (fNovfl)
1106     target.fNovfl = fNovfl;
1107   if (fData)
1108     target.fData = fData;
1109   if (fErr2)
1110     target.fErr2 = fErr2;
1111   
1112 }
1113 //____________________________________________________________________
1114 void AliCFGrid::SetExcludeOffEntriesInProj(Bool_t in)
1115 {
1116   //
1117   // require under/overflows in 'hidden dimensions' to be excluded
1118   // or included, when performing projections.
1119   // For AliCFGrid implementation, only option = kTRUE is available
1120
1121   if(!in){
1122     AliInfo(Form("This option is not available for AliCFGrid, Under/Overflows in hidden dimensions are always excluded")); 
1123     return;
1124   }
1125
1126   fExclOffEntriesInProj=in;
1127
1128 //____________________________________________________________________
1129 Bool_t AliCFGrid::GetExcludeOffEntriesInProj( ) const 
1130 {
1131   //
1132   // return flag saying whether under/overflows are excluded in projections 
1133   //
1134   
1135   return fExclOffEntriesInProj;
1136