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