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, Float_t min, Float_t max)
 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     Int_t binMin = GetMeasuredHistogram()->GetYaxis()->FindBin(min);
 116     Int_t binMax = GetMeasuredHistogram()->GetYaxis()->FindBin(max);
 117
 118     if (min==0 && max==0) {
 119       meas1D = GetMeasuredHistogram()->ProjectionX();
 120       gene1D = GetGeneratedHistogram()->ProjectionX();
 121     }
 122     else {
 123       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including y-bins %d to %d \n", binMin, binMax));
 124
 125       meas1D = GetMeasuredHistogram()->ProjectionX("pm",binMin,binMax);
 126       gene1D = GetGeneratedHistogram()->ProjectionX("pg",binMin,binMax);
 127     }
 128   }
 129   if (strcmp(opt,"y")==0) {
 130     Int_t binMin = GetMeasuredHistogram()->GetXaxis()->FindBin(min);
 131     Int_t binMax = GetMeasuredHistogram()->GetXaxis()->FindBin(max);
 132
 133     if (min==0 && max==0) {
 134       meas1D = GetMeasuredHistogram()->ProjectionY();
 135       gene1D = GetGeneratedHistogram()->ProjectionY();
 136     }
 137     else {
 138       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));
 139
 140       meas1D = GetMeasuredHistogram()->ProjectionY("pm", binMin, binMax);
 141       gene1D = GetGeneratedHistogram()->ProjectionY("pg", binMin, binMax);
 142     }
 143   }
 144   gene1D->Sumw2();
 145
 146   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
 147   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 148
 149   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
 150
 151   return (TH1F*)gene1D;
 152 }
 153
 154 //____________________________________________________________________
 155 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
 156 {
 157   // add value to measured histogram
 158   GetMeasuredHistogram()->Fill(ax, ay);
 159 }
 160
 161 //____________________________________________________________________
 162 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
 163 {
 164   // add value to generated histogram
 165   GetGeneratedHistogram()->Fill(ax, ay);
 166 }
 167
 168 //____________________________________________________________________
 169 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
 170 {
 171   // returns a value of the correction map
 172   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
 173 }
 174
 175 //____________________________________________________________________
 176 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 177 {
 178   // remove edges of correction histogram by removing
 179   // - bins with content less than cut
 180   // - bins next to bins with zero bin content
 181
 182   Int_t nBinsX = fhCorr->GetNbinsX();
 183   Int_t nBinsY = fhCorr->GetNbinsY();
 184
 185   // set bin content to zero for bins with content smaller cut
 186   for (Int_t bx=0; bx<=nBinsX; bx++) {
 187     for (Int_t by=0; by<=nBinsY; by++) {
 188       if (fhCorr->GetBinContent(bx,by)>cut) {
 189           fhCorr->SetBinContent(bx,by,0);
 190           fhCorr->SetBinError(bx,by,0);
 191       }
 192     }
 193   }
 194
 195   // set bin content to zero for bins next to bins with zero
 196   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 197   tmp->Reset();
 198
 199   Bool_t done = kFALSE;
 200   Int_t nBinsXCount = 0;
 201   Int_t nBinsYCount = 0;
 202   while (!done) {
 203     if (nBinsXCount<nBinsXedge)
 204       for (Int_t bx=0; bx<=nBinsX; bx++) {
 205         for (Int_t by=0; by<=nBinsY; by++) {
 206           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 207               (fhCorr->GetBinContent(bx-1,by)==0))
 208             tmp->SetBinContent(bx,by,1);
 209
 210         }
 211       }
 212     if (nBinsYCount<nBinsYedge)
 213       for (Int_t bx=0; bx<=nBinsX; bx++) {
 214         for (Int_t by=0; by<=nBinsY; by++) {
 215           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 216               (fhCorr->GetBinContent(bx,by-1)==0))
 217             tmp->SetBinContent(bx,by,1);
 218         }
 219       }
 220     for (Int_t bx=0; bx<=nBinsX; bx++) {
 221       for (Int_t by=0; by<=nBinsY; by++) {
 222         if (tmp->GetBinContent(bx,by)==1) {
 223           fhCorr->SetBinContent(bx,by,0);
 224           fhCorr->SetBinError(bx,by,0);
 225         }
 226       }
 227     }
 228     nBinsXCount++;
 229     nBinsYCount++;
 230     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 231   }
 232   tmp->Delete();
 233
 234 }
 235