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(const 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   TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
 167   divided->Reset();
 168
 169   divided->Divide(gene1D, meas1D, 1, 1, "B");
 170
 171   Printf("%p %p", gene1D, meas1D);
 172
 173   return (TH1F*)divided;
 174 }
 175
 176 //____________________________________________________________________
 177 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
 178 {
 179   // add value to measured histogram
 180   ((TH2F*)fhMeas)->Fill(ax, ay);
 181 }
 182
 183 //____________________________________________________________________
 184 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
 185 {
 186   // add value to generated histogram
 187   ((TH2F*)fhGene)->Fill(ax, ay);
 188 }
 189
 190 //____________________________________________________________________
 191 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
 192 {
 193   // returns a value of the correction map
 194   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
 195 }
 196
 197 //____________________________________________________________________
 198 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 199 {
 200   // remove edges of correction histogram by removing
 201   // - bins with content less than cut
 202   // - bins next to bins with zero bin content
 203
 204   Int_t nBinsX = fhCorr->GetNbinsX();
 205   Int_t nBinsY = fhCorr->GetNbinsY();
 206
 207   // set bin content to zero for bins with content smaller cut
 208   for (Int_t bx=0; bx<=nBinsX; bx++) {
 209     for (Int_t by=0; by<=nBinsY; by++) {
 210       if (fhCorr->GetBinContent(bx,by)>cut) {
 211           fhCorr->SetBinContent(bx,by,0);
 212           fhCorr->SetBinError(bx,by,0);
 213       }
 214     }
 215   }
 216
 217   // set bin content to zero for bins next to bins with zero
 218   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 219   tmp->Reset();
 220
 221   Bool_t done = kFALSE;
 222   Int_t nBinsXCount = 0;
 223   Int_t nBinsYCount = 0;
 224   while (!done) {
 225     if (nBinsXCount<nBinsXedge)
 226       for (Int_t bx=0; bx<=nBinsX; bx++) {
 227         for (Int_t by=0; by<=nBinsY; by++) {
 228           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 229               (fhCorr->GetBinContent(bx-1,by)==0))
 230             tmp->SetBinContent(bx,by,1);
 231
 232         }
 233       }
 234     if (nBinsYCount<nBinsYedge)
 235       for (Int_t bx=0; bx<=nBinsX; bx++) {
 236         for (Int_t by=0; by<=nBinsY; by++) {
 237           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 238               (fhCorr->GetBinContent(bx,by-1)==0))
 239             tmp->SetBinContent(bx,by,1);
 240         }
 241       }
 242     for (Int_t bx=0; bx<=nBinsX; bx++) {
 243       for (Int_t by=0; by<=nBinsY; by++) {
 244         if (tmp->GetBinContent(bx,by)==1) {
 245           fhCorr->SetBinContent(bx,by,0);
 246           fhCorr->SetBinError(bx,by,0);
 247         }
 248       }
 249     }
 250     nBinsXCount++;
 251     nBinsYCount++;
 252     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 253   }
 254   tmp->Delete();
 255
 256 }
 257
 258 //____________________________________________________________________
 259 void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
 260 {
 261         // rebins the histograms, recalculates the correction
 262
 263         GetGeneratedHistogram()->Rebin2D(x, y);
 264         GetMeasuredHistogram()->Rebin2D(x, y);
 265         GetCorrectionHistogram()->Rebin2D(x, y);
 266         Divide();
 267 }