]> git.uio.no Git - u/mrichter/AliRoot.git/blame - CORRFW/AliCFGrid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
further mods to include AliCFGridSparse, now included in compilation, and fixup for...
[u/mrichter/AliRoot.git] / CORRFW / AliCFGrid.cxx
CommitLineData
563113d0 1/* $Id$ */
1e9dad92 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16//---------------------------------------------------------------------//
17// //
18// AliCFGrid Class //
19// Class to accumulate data on an N-dimensional grid, to be used //
20// as input to get corrections for Reconstruction & Trigger efficiency //
21// The class uses a one-dimensional array of floats to store the grid //
22// --Author : S.Arcelli //
23// Still to be done: //
24// --Implement methods to merge cells //
25// --Interpolate among bins in a range //
26// This implementation will be aventually replaced byAliCFGridSparse //
27//---------------------------------------------------------------------//
563113d0 28//
29//
30#include <AliLog.h>
31#include "AliCFGrid.h"
32#include "TMath.h"
33#include "TROOT.h"
1e9dad92 34#include "TH1D.h"
35#include "TH2D.h"
36#include "TH3D.h"
563113d0 37
38//____________________________________________________________________
39ClassImp(AliCFGrid)
40
41//____________________________________________________________________
42AliCFGrid::AliCFGrid() :
1e9dad92 43 AliCFVGrid(),
563113d0 44 fNentriesTot(0),
45 fNunfl(0x0),
46 fNovfl(0x0),
47 fData(0x0),
48 fErr2(0x0)
49{
50 // default constructor
51}
52//____________________________________________________________________
53AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title) :
1e9dad92 54 AliCFVGrid(name,title),
563113d0 55 fNentriesTot(0),
56 fNunfl(0x0),
57 fNovfl(0x0),
58 fData(0x0),
59 fErr2(0x0)
60{
61 // default constructor
62}
63
64//____________________________________________________________________
1e9dad92 65AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title, const Int_t nVarIn, const Int_t * nBinIn, const Double_t *binLimitsIn) :
66 AliCFVGrid(name,title,nVarIn,nBinIn,binLimitsIn),
563113d0 67 fNentriesTot(0),
68 fNunfl(0x0),
69 fNovfl(0x0),
70 fData(0x0),
71 fErr2(0x0)
72{
73 //
74 // main constructor
75 //
1e9dad92 76
77 //The over/underflows
563113d0 78 fNunfl=new Float_t[fNVar];
79 fNovfl= new Float_t[fNVar];
80
81
563113d0 82 //Initialization
563113d0 83 fNentriesTot =0;
84 for(Int_t j=0;j<fNVar;j++){
85 fNunfl[j] =0;
86 fNovfl[j] =0;
87 }
1e9dad92 88
89
90 // the grid
91
92 fData = new Float_t[fNDim]; //num
93
94 //Initialization
95
563113d0 96}
97
98//____________________________________________________________________
99AliCFGrid::AliCFGrid(const AliCFGrid& c) :
1e9dad92 100 AliCFVGrid(c),
563113d0 101 fNentriesTot(0),
102 fNunfl(0x0),
103 fNovfl(0x0),
104 fData(0x0),
105 fErr2(0x0)
106{
107 //
108 // copy constructor
109 //
110 ((AliCFGrid &)c).Copy(*this);
111}
112
113//____________________________________________________________________
114AliCFGrid::~AliCFGrid()
115{
116 //
117 // destructor
118 //
1e9dad92 119 if(fNunfl)delete fNunfl;
120 if(fNovfl)delete fNovfl;
121 if(fData)delete fData;
563113d0 122 if(fSumW2)delete fErr2;
563113d0 123
124}
563113d0 125//____________________________________________________________________
126AliCFGrid &AliCFGrid::operator=(const AliCFGrid &c)
127{
128 //
129 // assigment operator
130 //
131 if (this != &c)
132 ((AliCFGrid &) c).Copy(*this);
133 return *this;
134}
135//____________________________________________________________________
136Float_t AliCFGrid::GetElement(Int_t bin) const
137{
138 //
139 // Returns grid element bin
140 //
141 if(bin>=fNDim){
142 AliInfo(Form(" element index outside the grid, return -1"));
143 return -1.;
144 }
145 return fData[bin];
146}
147//____________________________________________________________________
148Float_t AliCFGrid::GetElement(Int_t *bin) const
149{
150 //
151 // Get the content in a bin corresponding to a set of bin indexes
152 //
1e9dad92 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);
563113d0 157}
158//____________________________________________________________________
1e9dad92 159Float_t AliCFGrid::GetElement(Double_t *var) const
563113d0 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
1e9dad92 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
563113d0 193 if(!(ovfl==1 || unfl==1)){
194 return GetElement(fIndex);
195 }
196 else{
197 AliInfo(Form(" input variables outside the grid, return -1"));
198 return -1.;
199 }
200}
201//____________________________________________________________________
202Float_t AliCFGrid::GetElementError(Int_t iel) const
203{
204 //
1e9dad92 205 // Return the error on grid element iel
563113d0 206 //
207 if(iel>=fNDim){
208 AliInfo(Form(" element index outside the grid, return -1"));
209 return -1.;
210 }
1e9dad92 211 if(fSumW2)return TMath::Sqrt(fErr2[iel]);
212 return TMath::Sqrt(fData[iel]);
563113d0 213}
214//____________________________________________________________________
215Float_t AliCFGrid::GetElementError(Int_t *bin) const
216{
217 //
1e9dad92 218 // Get the error in a bin corresponding to a set of bin indeces
563113d0 219 //
1e9dad92 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);
563113d0 224
225}
226//____________________________________________________________________
1e9dad92 227Float_t AliCFGrid::GetElementError(Double_t *var) const
563113d0 228{
229 //
1e9dad92 230 // Get the error in a bin corresponding to a set of input variables
563113d0 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
1e9dad92 256 //move to the TH/THnSparse convention in N-dim bin numbering
257 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
563113d0 258 if(!(ovfl==1 || unfl==1)){
259 return GetElementError(fIndex);
260 }
261 else{
262 AliInfo(Form(" input variables outside the grid, return -1"));
263 return -1.;
264 }
265}
266//____________________________________________________________________
267void AliCFGrid::SetElement(Int_t iel, Float_t val)
268{
269 //
270 // Sets grid element iel to val
271 //
272 if(iel>=fNDim){
273 AliInfo(Form(" element index outside the grid, no value set"));
274 }else {
275 fData[iel]=val;
276 }
277}
278//____________________________________________________________________
279void AliCFGrid::SetElement(Int_t *bin, Float_t val)
280{
281 //
282 // Sets grid element of bin indeces bin to val
283 //
1e9dad92 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);
563113d0 288}
289//____________________________________________________________________
1e9dad92 290void AliCFGrid::SetElement(Double_t *var, Float_t val)
563113d0 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
1e9dad92 319 //move to the TH/THnSparse convention in N-dim bin numbering
320 for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
563113d0 321 if(!(ovfl==1 || unfl==1)){
322 SetElement(fIndex,val);
323 }
324 else{
325 AliInfo(Form(" input variables outside the grid, no value set"));
326 }
327}
328//____________________________________________________________________
329void AliCFGrid::SetElementError(Int_t iel, Float_t val)
330{
331 //
1e9dad92 332 // Set squared error on grid element iel to val*val
563113d0 333 //
334 if(iel>=fNDim){
335 AliInfo(Form(" element index outside the grid, no value set"));
336 return;
337 }
338 if(!fErr2)SumW2();
1e9dad92 339 fErr2[iel]=val*val;
563113d0 340}
341//____________________________________________________________________
342void AliCFGrid::SetElementError(Int_t *bin, Float_t val)
343{
344 //
345 // Set squared error to val on grid element of bin indeces bin
346 //
1e9dad92 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);
563113d0 351}
352//____________________________________________________________________
1e9dad92 353void AliCFGrid::SetElementError(Double_t *var, Float_t val)
563113d0 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
1e9dad92 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
563113d0 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//____________________________________________________________________
1e9dad92 393void AliCFGrid::Fill(Double_t *var, Double_t weight)
563113d0 394{
1e9dad92 395
563113d0 396 //
1e9dad92 397 // Fill the grid,
398 // given a set of values of the input variable,
399 // with weight (by default w=1)
563113d0 400 //
1e9dad92 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;
563113d0 411 }
1e9dad92 412
563113d0 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]){
1e9dad92 424 unfl[i]=1;
425 isunfl=1;
563113d0 426 }
427
428 //overflows
429
430 if(var[i] > bins[nbins-1]){
1e9dad92 431 ovfl[i]=1;
432 isovfl=1;
563113d0 433 }
434 delete [] bins;
435 }
436
1e9dad92 437 //exclusive under/overflows
563113d0 438
563113d0 439 for(Int_t i=0;i<fNVar;i++){
1e9dad92 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;
563113d0 444 }
1e9dad92 445 if(add && unfl[i]==1)fNunfl[i]++;
446 if(add && ovfl[i]==1)fNovfl[i]++;
563113d0 447 }
448
1e9dad92 449 delete [] unfl;
450 delete [] ovfl;
451
563113d0 452 // Total number of entries, overflows and underflows
453
454 fNentriesTot++;
563113d0 455
456 //if not ovfl/unfl, fill the element
1e9dad92 457
458 if(!(isovfl==1 || isunfl==1)){
563113d0 459 Int_t ind =GetBinIndex(fIndex);
460 fData[ind]+=weight;
461 if(fSumW2)fErr2[ind]+=(weight*weight);
462 }
463}
464//____________________________________________________________________
465Float_t AliCFGrid::GetOverFlows( Int_t ivar) const {
466 //
467 // Get overflows in variable var
468 //
469 return fNovfl[ivar];
470}
471//____________________________________________________________________
563113d0 472Float_t AliCFGrid::GetUnderFlows( Int_t ivar) const {
473 //
474 // Get overflows in variable var
475 //
476 return fNunfl[ivar];
477}
478//____________________________________________________________________
563113d0 479Float_t AliCFGrid::GetEntries( ) const {
480 //
481 // Get total entries (in grid + over/underflows)
482 //
483 return fNentriesTot;
484}
563113d0 485//___________________________________________________________________
1e9dad92 486TH1D *AliCFGrid::Project(Int_t ivar) const
563113d0 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
1e9dad92 505 TH1D *proj1D=0;
563113d0 506
507 //check if a projection with identical name exist
508 TObject *obj = gROOT->FindObject(pname);
1e9dad92 509 if (obj && obj->InheritsFrom("TH1D")) {
510 proj1D = (TH1D*)obj;
563113d0 511 proj1D->Reset();
512 }
513
514 if(!proj1D){
1e9dad92 515 proj1D =new TH1D(pname,htitle, nbins, bins);
563113d0 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));
1e9dad92 541 proj1D->SetEntries(sum+GetUnderFlows(ivar)+GetOverFlows(ivar));
563113d0 542 return proj1D;
543}
544
545//___________________________________________________________________
1e9dad92 546TH2D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2) const
563113d0 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
1e9dad92 569 TH2D *proj2D=0;
563113d0 570
571 //check if a projection with identical name exist
572 TObject *obj = gROOT->FindObject(pname);
1e9dad92 573 if (obj && obj->InheritsFrom("TH2D")) {
574 proj2D = (TH2D*)obj;
563113d0 575 proj2D->Reset();
576 }
577
578 if(!proj2D){
1e9dad92 579 proj2D =new TH2D(pname,htitle, nbins1, bins1,nbins2,bins2);
563113d0 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));
1e9dad92 623 proj2D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2));
563113d0 624 return proj2D;
625}
626//___________________________________________________________________
1e9dad92 627TH3D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
563113d0 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
1e9dad92 655 TH3D *proj3D=0;
563113d0 656
657 //check if a projection with identical name exist
658 TObject *obj = gROOT->FindObject(pname);
1e9dad92 659 if (obj && obj->InheritsFrom("TH3D")) {
660 proj3D = (TH3D*)obj;
563113d0 661 proj3D->Reset();
662 }
663
664 if(!proj3D){
1e9dad92 665 proj3D =new TH3D(pname,htitle, nbins1,bins1,nbins2,bins2,nbins3,bins3);
563113d0 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;
1e9dad92 719 proj3D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2)+GetUnderFlows(ivar3)+GetOverFlows(ivar3));
563113d0 720 return proj3D;
721}
722
723//___________________________________________________________________
1e9dad92 724TH1D *AliCFGrid::Slice(Int_t ivar, Double_t *varMin, Double_t* varMax) const
563113d0 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
1e9dad92 743 TH1D *proj1D=0;
563113d0 744
745 //check if a projection with identical name exist
746 TObject *obj = gROOT->FindObject(pname);
1e9dad92 747 if (obj && obj->InheritsFrom("TH1D")) {
748 proj1D = (TH1D*)obj;
563113d0 749 proj1D->Reset();
750 }
751
752 if(!proj1D){
1e9dad92 753 proj1D =new TH1D(pname,htitle, nbins, bins);
563113d0 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);
563113d0 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
563113d0 810
811//____________________________________________________________________
1e9dad92 812void AliCFGrid::Add(AliCFVGrid* aGrid, Double_t c)
563113d0 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){
1e9dad92 832 Float_t err=aGrid->GetElementError(iel);
563113d0 833 fErr2[iel]+=c*c*err*err;
834 }
835 }
836
837 //Add entries, overflows and underflows
838
839 fNentriesTot+= c*aGrid->GetEntries();
563113d0 840 for(Int_t j=0;j<fNVar;j++){
841 fNunfl[j]+= c*aGrid->GetUnderFlows(j);
db6722a5 842 fNovfl[j]+= c*aGrid->GetOverFlows(j);
563113d0 843 }
844}
845//____________________________________________________________________
1e9dad92 846void AliCFGrid::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
563113d0 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){
1e9dad92 870 err1=aGrid1->GetElementError(iel);
871 err2=aGrid2->GetElementError(iel);
872 SetElementError(iel,TMath::Sqrt(c1*c1*err1*err1+c2*c2*err2*err2));
563113d0 873 }
874 }
875
876 //Add entries, overflows and underflows
877
878 fNentriesTot= c1*aGrid1->GetEntries()+c2*aGrid2->GetEntries();
563113d0 879 for(Int_t j=0;j<fNVar;j++){
880 fNunfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
db6722a5 881 fNovfl[j]= c1*aGrid1->GetOverFlows(j)+c2*aGrid2->GetOverFlows(j);
563113d0 882 }
883}
884//____________________________________________________________________
1e9dad92 885void AliCFGrid::Multiply(AliCFVGrid* aGrid, Double_t c)
563113d0 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){
1e9dad92 909 err1=GetElementError(iel);
910 err2=aGrid->GetElementError(iel);
911 SetElementError(iel,TMath::Sqrt(c*c*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
563113d0 912 }
913 }
914
915 //Set entries to the number of bins, preserve original overflows and underflows
916
917 fNentriesTot=fNDim;
563113d0 918 for(Int_t j=0;j<fNVar;j++){
919 fNunfl[j]= GetUnderFlows(j);
db6722a5 920 fNovfl[j]= GetOverFlows(j);
563113d0 921 }
922}
923//____________________________________________________________________
1e9dad92 924void AliCFGrid::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1, Double_t c2)
563113d0 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){
1e9dad92 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)));
563113d0 950 }
951 }
952
953 //Set entries to the number of bins, preserve original overflows and underflows
954
955 fNentriesTot=fNDim;
563113d0 956 for(Int_t j=0;j<fNVar;j++){
957 fNunfl[j]= GetUnderFlows(j);
db6722a5 958 fNovfl[j]= GetOverFlows(j);
563113d0 959 }
960}
961//____________________________________________________________________
1e9dad92 962void AliCFGrid::Divide(AliCFVGrid* aGrid, Double_t c)
563113d0 963{
964 //
965 //divide current grid by grid aGrid
966 //
967
563113d0 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
1e9dad92 980 Float_t cont1,cont2,err1,err2,den;
563113d0 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){
1e9dad92 987 err1=GetElementError(iel);
988 err2=aGrid->GetElementError(iel);
563113d0 989 if(!cont2){SetElementError(iel,0.); continue;}
1e9dad92 990 den=cont2*cont2*cont2*c*c;
991 SetElementError(iel,TMath::Sqrt((cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
563113d0 992 }
993 }
1e9dad92 994
563113d0 995 //Set entries to the number of bins, preserve original overflows and underflows
996
997 fNentriesTot=fNDim;
563113d0 998 for(Int_t j=0;j<fNVar;j++){
999 fNunfl[j]= GetUnderFlows(j);
db6722a5 1000 fNovfl[j]= GetOverFlows(j);
563113d0 1001 }
1002}
1003//____________________________________________________________________
1e9dad92 1004void AliCFGrid::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
563113d0 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){
1e9dad92 1033 err1=aGrid1->GetElementError(iel);
1034 err2=aGrid2->GetElementError(iel);
563113d0 1035 if(!cont2){SetElementError(iel,0.); continue;}
1036 if (opt.Contains("B")){
1037 if(cont1!=cont2){
1038 r=cont1/cont2;
1e9dad92 1039 SetElementError(iel,TMath::Sqrt(TMath::Abs(((1.-2.*r)*err1*err1+r*r*err2*err2)/(cont2*cont2))));
563113d0 1040 }else{
1041 SetElementError(iel,0.);
1042 }
1043 }else{
1044 den=cont2*cont2*cont2*cont2*c2*c2;
1e9dad92 1045 SetElementError(iel,TMath::Sqrt(c1*c1*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
563113d0 1046 }
1047 }
1048 }
1049
1050 //Set entries to the number of bins, preserve original overflows and underflows
1051
1052 fNentriesTot=fNDim;
563113d0 1053 for(Int_t j=0;j<fNVar;j++){
1054 fNunfl[j]= GetUnderFlows(j);
db6722a5 1055 fNovfl[j]= GetOverFlows(j);
563113d0 1056 }
1057}
1058//____________________________________________________________________
1059void AliCFGrid::SumW2()
1060{
1061 //
1062 //set calculation of the squared sum of the weighted entries
1063 //
1064 if(!fSumW2){
1065 fErr2=new Float_t [fNDim];
1066 //init....
1067 for(Int_t iel=0;iel<fNDim;iel++){
1068 fErr2[iel]=fData[iel];
1069 }
1070 }
1071
1072 fSumW2=kTRUE;
1073}
563113d0 1074//____________________________________________________________________
1075void AliCFGrid::Copy(TObject& c) const
1076{
1077 //
1078 // copy function
1079 //
1080 AliCFGrid& target = (AliCFGrid &) c;
1081
563113d0 1082 target.fNentriesTot = fNentriesTot;
563113d0 1083 if (fNunfl)
1084 target.fNunfl = fNunfl;
1085 if (fNovfl)
1086 target.fNovfl = fNovfl;
563113d0 1087 if (fData)
1088 target.fData = fData;
1089 if (fErr2)
1090 target.fErr2 = fErr2;
1091
1092}