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",title),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  60   fhGene  = new TH2F("generated",  Form("%s generated",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  61   fhCorr  = new TH2F("correction", Form("%s correction",title), 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 TH2F* AliCorrectionMatrix2D::GetGeneratedHistogram() const
 105 {
 106   // return generated histogram casted to correct type
 107   return dynamic_cast<TH2F*> (fhGene);
 108 }
 109
 110 TH2F* AliCorrectionMatrix2D::GetMeasuredHistogram() const
 111 {
 112   // return measured histogram casted to correct type
 113   return dynamic_cast<TH2F*> (fhMeas);
 114 }
 115
 116 //____________________________________________________________________
 117 TH1F* 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       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including y-bins %d to %d \n", binMin, binMax));
 136
 137       meas1D = ((TH2F*)fhMeas)->ProjectionX("pm",binMin,binMax);
 138       gene1D = ((TH2F*)fhGene)->ProjectionX("pg",binMin,binMax);
 139     }
 140   }
 141   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       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));
 151
 152       meas1D = ((TH2F*)fhMeas)->ProjectionY("pm", binMin, binMax);
 153       gene1D = ((TH2F*)fhGene)->ProjectionY("pg", binMin, binMax);
 154     }
 155   }
 156   gene1D->Sumw2();
 157
 158   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
 159   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 160
 161   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
 162
 163   return (TH1F*)gene1D;
 164 }
 165
 166 //____________________________________________________________________
 167 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
 168 {
 169   // add value to measured histogram
 170   ((TH2F*)fhMeas)->Fill(ax, ay);
 171 }
 172
 173 //____________________________________________________________________
 174 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
 175 {
 176   // add value to generated histogram
 177   ((TH2F*)fhGene)->Fill(ax, ay);
 178 }
 179
 180 //____________________________________________________________________
 181 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
 182 {
 183   // returns a value of the correction map
 184   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
 185 }
 186
 187 //____________________________________________________________________
 188 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 189 {
 190   // remove edges of correction histogram by removing
 191   // - bins with content less than cut
 192   // - bins next to bins with zero bin content
 193
 194   Int_t nBinsX = fhCorr->GetNbinsX();
 195   Int_t nBinsY = fhCorr->GetNbinsY();
 196
 197   // set bin content to zero for bins with content smaller cut
 198   for (Int_t bx=0; bx<=nBinsX; bx++) {
 199     for (Int_t by=0; by<=nBinsY; by++) {
 200       if (fhCorr->GetBinContent(bx,by)>cut) {
 201           fhCorr->SetBinContent(bx,by,0);
 202           fhCorr->SetBinError(bx,by,0);
 203       }
 204     }
 205   }
 206
 207   // set bin content to zero for bins next to bins with zero
 208   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 209   tmp->Reset();
 210
 211   Bool_t done = kFALSE;
 212   Int_t nBinsXCount = 0;
 213   Int_t nBinsYCount = 0;
 214   while (!done) {
 215     if (nBinsXCount<nBinsXedge)
 216       for (Int_t bx=0; bx<=nBinsX; bx++) {
 217         for (Int_t by=0; by<=nBinsY; by++) {
 218           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 219               (fhCorr->GetBinContent(bx-1,by)==0))
 220             tmp->SetBinContent(bx,by,1);
 221
 222         }
 223       }
 224     if (nBinsYCount<nBinsYedge)
 225       for (Int_t bx=0; bx<=nBinsX; bx++) {
 226         for (Int_t by=0; by<=nBinsY; by++) {
 227           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 228               (fhCorr->GetBinContent(bx,by-1)==0))
 229             tmp->SetBinContent(bx,by,1);
 230         }
 231       }
 232     for (Int_t bx=0; bx<=nBinsX; bx++) {
 233       for (Int_t by=0; by<=nBinsY; by++) {
 234         if (tmp->GetBinContent(bx,by)==1) {
 235           fhCorr->SetBinContent(bx,by,0);
 236           fhCorr->SetBinError(bx,by,0);
 237         }
 238       }
 239     }
 240     nBinsXCount++;
 241     nBinsYCount++;
 242     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 243   }
 244   tmp->Delete();
 245
 246 }
 247