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(const Char_t* name, const Char_t* title,
  36                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
  37                                        Int_t nBinY, Float_t Ymin, Float_t Ymax)
  38   : AliCorrectionMatrix(name, title)
  39 {
  40   //
  41   // constructor
  42   //
  43
  44   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  45   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  46   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  47
  48   fhMeas->Sumw2();
  49   fhGene->Sumw2();
  50   fhCorr->Sumw2();
  51 }
  52
  53 //____________________________________________________________________
  54 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
  55                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y)
  56   : AliCorrectionMatrix(name, title)
  57 {
  58   //
  59   // constructor
  60   //
  61
  62   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
  63   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
  64   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
  65
  66   fhMeas->Sumw2();
  67   fhGene->Sumw2();
  68   fhCorr->Sumw2();
  69 }
  70
  71 //____________________________________________________________________
  72 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
  73 {
  74   //
  75   // destructor
  76   //
  77
  78   // histograms already deleted in base class
  79 }
  80
  81 TH2F* AliCorrectionMatrix2D::GetGeneratedHistogram() const
  82 {
  83   // return generated histogram casted to correct type
  84   return dynamic_cast<TH2F*> (fhGene);
  85 }
  86
  87 TH2F* AliCorrectionMatrix2D::GetMeasuredHistogram() const
  88 {
  89   // return measured histogram casted to correct type
  90   return dynamic_cast<TH2F*> (fhMeas);
  91 }
  92
  93 //____________________________________________________________________
  94 TH1F* AliCorrectionMatrix2D::Get1DCorrection(Char_t* opt)
  95 {
  96   //
  97   // integrate the correction over one variable
  98   //
  99
 100   TH1D* meas1D = 0;
 101   TH1D* gene1D = 0;
 102
 103   if (strcmp(opt,"x")==0) {
 104     meas1D = GetMeasuredHistogram()->ProjectionX();
 105     gene1D = GetGeneratedHistogram()->ProjectionX();
 106   }
 107   if (strcmp(opt,"y")==0) {
 108     meas1D = GetMeasuredHistogram()->ProjectionY();
 109     gene1D = GetGeneratedHistogram()->ProjectionY();
 110   }
 111   gene1D->Sumw2();
 112
 113   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
 114   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 115
 116   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
 117
 118   return (TH1F*)gene1D;
 119 }
 120
 121 //____________________________________________________________________
 122 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
 123 {
 124   // add value to measured histogram
 125   GetMeasuredHistogram()->Fill(ax, ay);
 126 }
 127
 128 //____________________________________________________________________
 129 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
 130 {
 131   // add value to generated histogram
 132   GetGeneratedHistogram()->Fill(ax, ay);
 133 }
 134
 135 //____________________________________________________________________
 136 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
 137 {
 138   // returns a value of the correction map
 139   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
 140 }
 141
 142 //____________________________________________________________________
 143 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 144 {
 145   // remove edges of correction histogram by removing
 146   // - bins with content less than cut
 147   // - bins next to bins with zero bin content
 148
 149   Int_t nBinsX = fhCorr->GetNbinsX();
 150   Int_t nBinsY = fhCorr->GetNbinsY();
 151
 152   // set bin content to zero for bins with content smaller cut
 153   for (Int_t bx=0; bx<=nBinsX; bx++) {
 154     for (Int_t by=0; by<=nBinsY; by++) {
 155       if (fhCorr->GetBinContent(bx,by)>cut) {
 156           fhCorr->SetBinContent(bx,by,0);
 157           fhCorr->SetBinError(bx,by,0);
 158       }
 159     }
 160   }
 161
 162   // set bin content to zero for bins next to bins with zero
 163   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 164   tmp->Reset();
 165
 166   Bool_t done = kFALSE;
 167   Int_t nBinsXCount = 0;
 168   Int_t nBinsYCount = 0;
 169   while (!done) {
 170     if (nBinsXCount<nBinsXedge)
 171       for (Int_t bx=0; bx<=nBinsX; bx++) {
 172         for (Int_t by=0; by<=nBinsY; by++) {
 173           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 174               (fhCorr->GetBinContent(bx-1,by)==0))
 175             tmp->SetBinContent(bx,by,1);
 176
 177         }
 178       }
 179     if (nBinsYCount<nBinsYedge)
 180       for (Int_t bx=0; bx<=nBinsX; bx++) {
 181         for (Int_t by=0; by<=nBinsY; by++) {
 182           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 183               (fhCorr->GetBinContent(bx,by-1)==0))
 184             tmp->SetBinContent(bx,by,1);
 185         }
 186       }
 187     for (Int_t bx=0; bx<=nBinsX; bx++) {
 188       for (Int_t by=0; by<=nBinsY; by++) {
 189         if (tmp->GetBinContent(bx,by)==1) {
 190           fhCorr->SetBinContent(bx,by,0);
 191           fhCorr->SetBinError(bx,by,0);
 192         }
 193       }
 194     }
 195     nBinsXCount++;
 196     nBinsYCount++;
 197     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 198   }
 199   tmp->Delete();
 200
 201 }
 202