PWG0/AliCorrectionMatrix2D.cxx

   1 /* $Id$ */
   2
   3 // ------------------------------------------------------
   4 //
   5 // Class to handle 2d-corrections.
   6 //
   7 // ------------------------------------------------------
   8 //
   9
  10 #include <TH2F.h>
  11 #include <TMath.h>
  12
  13 #include <AliLog.h>
  14
  15 #include "AliCorrectionMatrix2D.h"
  16
  17 //____________________________________________________________________
  18 ClassImp(AliCorrectionMatrix2D)
  19
  20 //____________________________________________________________________
  21 AliCorrectionMatrix2D::AliCorrectionMatrix2D() :
  22   AliCorrectionMatrix()
  23 {
  24   // default constructor
  25 }
  26
  27 //____________________________________________________________________
  28 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const AliCorrectionMatrix2D& c)
  29   : AliCorrectionMatrix(c)
  30 {
  31   // copy constructor
  32   ((AliCorrectionMatrix2D &)c).Copy(*this);
  33 }
  34
  35 //____________________________________________________________________
  36 AliCorrectionMatrix2D &AliCorrectionMatrix2D::operator=(const AliCorrectionMatrix2D &c)
  37 {
  38   // assigment operator
  39
  40   if (this != &c)
  41     ((AliCorrectionMatrix2D &) c).Copy(*this);
  42
  43   return *this;
  44 }
  45
  46 //____________________________________________________________________
  47 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
  48                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
  49                                        Int_t nBinY, Float_t Ymin, Float_t Ymax)
  50   : AliCorrectionMatrix(name, title)
  51 {
  52   //
  53   // constructor
  54   //
  55
  56   // do not add this hists to the directory
  57   Bool_t oldStatus = TH1::AddDirectoryStatus();
  58   TH1::AddDirectory(kFALSE);
  59
  60   fhMeas  = new TH2F("measured",   Form("%s measured", GetTitle()),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  61   fhGene  = new TH2F("generated",  Form("%s generated", GetTitle()),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  62   fhCorr  = new TH2F("correction", Form("%s correction", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  63
  64   TH1::AddDirectory(oldStatus);
  65
  66   fhMeas->Sumw2();
  67   fhGene->Sumw2();
  68   fhCorr->Sumw2();
  69 }
  70
  71 //____________________________________________________________________
  72 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
  73                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y)
  74   : AliCorrectionMatrix(name, title)
  75 {
  76   //
  77   // constructor
  78   //
  79
  80   // do not add this hists to the directory
  81   Bool_t oldStatus = TH1::AddDirectoryStatus();
  82   TH1::AddDirectory(kFALSE);
  83
  84         fhMeas  = new TH2F("measured",   Form("%s measured",title),   nBinX, X, nBinY, Y);
  85   fhGene  = new TH2F("generated",  Form("%s generated",title),  nBinX, X, nBinY, Y);
  86   fhCorr  = new TH2F("correction", Form("%s correction",title), nBinX, X, nBinY, Y);
  87
  88   TH1::AddDirectory(oldStatus);
  89
  90         fhMeas->Sumw2();
  91   fhGene->Sumw2();
  92   fhCorr->Sumw2();
  93 }
  94
  95 //____________________________________________________________________
  96 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
  97 {
  98   //
  99   // destructor
 100   //
 101
 102   // histograms already deleted in base class
 103 }
 104
 105 TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
 106 {
 107   // return generated histogram casted to correct type
 108   return dynamic_cast<TH2*> (fhGene);
 109 }
 110
 111 TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
 112 {
 113   // return measured histogram casted to correct type
 114   return dynamic_cast<TH2*> (fhMeas);
 115 }
 116
 117 //____________________________________________________________________
 118 TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const Char_t* opt, Float_t min, Float_t max, Bool_t binomialErrors)
 119 {
 120   //
 121   // integrate the correction over one variable
 122   //
 123
 124   TH1D* meas1D = 0;
 125   TH1D* gene1D = 0;
 126
 127   if (strcmp(opt,"x")==0) {
 128     Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
 129     Int_t binMax = fhGene->GetYaxis()->FindBin(max);
 130
 131     if (min>=max) {
 132       meas1D = ((TH2F*)fhMeas)->ProjectionX();
 133       gene1D = ((TH2F*)fhGene)->ProjectionX();
 134     }
 135     else {
 136       Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
 137
 138       meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
 139       gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
 140     }
 141   }
 142   else if (strcmp(opt,"y")==0) {
 143     Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
 144     Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
 145
 146     if (min>=max) {
 147       meas1D = ((TH2F*)fhMeas)->ProjectionY();
 148       gene1D = ((TH2F*)fhGene)->ProjectionY();
 149     }
 150     else {
 151       Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);
 152
 153       meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
 154       gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
 155     }
 156   }
 157   else {
 158     Printf("ERROR: Invalid option");
 159     return 0;
 160   }
 161
 162   if (!binomialErrors)
 163   {
 164     // set the errors on gene manually, and clear the ones on meas.
 165     gene1D->Sumw2();
 166     for (Int_t bin=0; bin <= gene1D->GetNbinsX()+1; bin++)
 167     {
 168       gene1D->SetBinError(bin, TMath::Sqrt(gene1D->GetBinContent(bin)));
 169       meas1D->SetBinError(bin, 0);
 170     }
 171   }
 172
 173   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
 174   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 175
 176   TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
 177   divided->Reset();
 178
 179   divided->Divide(gene1D, meas1D, 1, 1, (binomialErrors) ? "B" : "");
 180
 181   Printf("%p %p", gene1D, meas1D);
 182
 183   return (TH1F*)divided;
 184 }
 185
 186 //____________________________________________________________________
 187 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
 188 {
 189   // add value to measured histogram
 190   ((TH2F*)fhMeas)->Fill(ax, ay);
 191 }
 192
 193 //____________________________________________________________________
 194 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
 195 {
 196   // add value to generated histogram
 197   ((TH2F*)fhGene)->Fill(ax, ay);
 198 }
 199
 200 //____________________________________________________________________
 201 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
 202 {
 203   // returns a value of the correction map
 204   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
 205 }
 206
 207 //____________________________________________________________________
 208 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 209 {
 210   // remove edges of correction histogram by removing
 211   // - bins with content less than cut
 212   // - bins next to bins with zero bin content
 213
 214   Int_t nBinsX = fhCorr->GetNbinsX();
 215   Int_t nBinsY = fhCorr->GetNbinsY();
 216
 217   // set bin content to zero for bins with content smaller cut
 218   for (Int_t bx=0; bx<=nBinsX; bx++) {
 219     for (Int_t by=0; by<=nBinsY; by++) {
 220       if (fhCorr->GetBinContent(bx,by)>cut) {
 221           fhCorr->SetBinContent(bx,by,0);
 222           fhCorr->SetBinError(bx,by,0);
 223       }
 224     }
 225   }
 226
 227   // set bin content to zero for bins next to bins with zero
 228   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 229   tmp->Reset();
 230
 231   Bool_t done = kFALSE;
 232   Int_t nBinsXCount = 0;
 233   Int_t nBinsYCount = 0;
 234   while (!done) {
 235     if (nBinsXCount<nBinsXedge)
 236       for (Int_t bx=0; bx<=nBinsX; bx++) {
 237         for (Int_t by=0; by<=nBinsY; by++) {
 238           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 239               (fhCorr->GetBinContent(bx-1,by)==0))
 240             tmp->SetBinContent(bx,by,1);
 241
 242         }
 243       }
 244     if (nBinsYCount<nBinsYedge)
 245       for (Int_t bx=0; bx<=nBinsX; bx++) {
 246         for (Int_t by=0; by<=nBinsY; by++) {
 247           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 248               (fhCorr->GetBinContent(bx,by-1)==0))
 249             tmp->SetBinContent(bx,by,1);
 250         }
 251       }
 252     for (Int_t bx=0; bx<=nBinsX; bx++) {
 253       for (Int_t by=0; by<=nBinsY; by++) {
 254         if (tmp->GetBinContent(bx,by)==1) {
 255           fhCorr->SetBinContent(bx,by,0);
 256           fhCorr->SetBinError(bx,by,0);
 257         }
 258       }
 259     }
 260     nBinsXCount++;
 261     nBinsYCount++;
 262     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 263   }
 264   tmp->Delete();
 265
 266 }
 267
 268 //____________________________________________________________________
 269 void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
 270 {
 271         // rebins the histograms, recalculates the correction
 272
 273         GetGeneratedHistogram()->Rebin2D(x, y);
 274         GetMeasuredHistogram()->Rebin2D(x, y);
 275         GetCorrectionHistogram()->Rebin2D(x, y);
 276         Divide();
 277 }