17f16dca86c241520b84681b8b8d0468c3369450
[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]; //num
93
94   //Initialization
95  
96 }
97
98 //____________________________________________________________________
99 AliCFGrid::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 //____________________________________________________________________
114 AliCFGrid::~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 //____________________________________________________________________
126 AliCFGrid &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 //____________________________________________________________________
136 Float_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 //____________________________________________________________________
148 Float_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 //____________________________________________________________________
159 Float_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 //____________________________________________________________________
202 Float_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 //____________________________________________________________________
215 Float_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 //____________________________________________________________________
227 Float_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 //____________________________________________________________________
267 void 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 //____________________________________________________________________
279 void 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 //____________________________________________________________________
290 void 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 //____________________________________________________________________
329 void 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 //____________________________________________________________________
342 void 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 //____________________________________________________________________
353 void 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 //____________________________________________________________________
393 void 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 //____________________________________________________________________
465 Float_t AliCFGrid::GetOverFlows( Int_t ivar) const {
466   //
467   // Get overflows in variable var 
468   //
469   return fNovfl[ivar];
470
471 //____________________________________________________________________
472 Float_t AliCFGrid::GetUnderFlows( Int_t ivar) const {
473   //
474   // Get overflows in variable var 
475   //
476   return fNunfl[ivar];
477
478 //____________________________________________________________________
479 Float_t AliCFGrid::GetEntries( ) const {
480   //
481   // Get total entries (in grid + over/underflows) 
482   //
483   return fNentriesTot;
484
485 //___________________________________________________________________
486 TH1D *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 //___________________________________________________________________
546 TH2D *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 //___________________________________________________________________
627 TH3D *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 //___________________________________________________________________
724 TH1D *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 //____________________________________________________________________
812 void 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 //____________________________________________________________________
846 void 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 //____________________________________________________________________
885 void 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 //____________________________________________________________________
924 void 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 //____________________________________________________________________
962 void 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 //____________________________________________________________________
1004 void 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 //____________________________________________________________________
1059 void 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 //____________________________________________________________________
1075 void 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 }