PWG0/AliCorrectionMatrix2D.cxx

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