PWG0/AliCorrectionMatrix3D.cxx

   1 /* $Id$ */
   2
   3 // ------------------------------------------------------
   4 //
   5 // Class to handle 3d-corrections.
   6 //
   7 // ------------------------------------------------------
   8 //
   9
  10 #include <TH3F.h>
  11 #include <TH2F.h>
  12 #include <TH1F.h>
  13 #include <TString.h>
  14
  15 #include <AliLog.h>
  16
  17 #include "AliCorrectionMatrix3D.h"
  18 #include "AliCorrectionMatrix2D.h"
  19 #include "AliPWG0Helper.h"
  20
  21 //____________________________________________________________________
  22 ClassImp(AliCorrectionMatrix3D)
  23
  24 //____________________________________________________________________
  25 AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
  26   AliCorrectionMatrix()
  27 {
  28   // default constructor
  29 }
  30
  31 //____________________________________________________________________
  32 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
  33   : AliCorrectionMatrix(c)
  34 {
  35   // copy constructor
  36   ((AliCorrectionMatrix3D &)c).Copy(*this);
  37 }
  38
  39 //____________________________________________________________________
  40 AliCorrectionMatrix3D &AliCorrectionMatrix3D::operator=(const AliCorrectionMatrix3D &c)
  41 {
  42   // assigment operator
  43
  44   if (this != &c)
  45     ((AliCorrectionMatrix3D &) c).Copy(*this);
  46
  47   return *this;
  48 }
  49
  50 //____________________________________________________________________
  51 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
  52               Int_t nBinX, Float_t Xmin, Float_t Xmax,
  53               Int_t nBinY, Float_t Ymin, Float_t Ymax,
  54               Int_t nBinZ, Float_t Zmin, Float_t Zmax)
  55   : AliCorrectionMatrix(name, title)
  56 {
  57   //
  58   // constructor
  59   //
  60
  61   Float_t* binLimitsX = new Float_t[nBinX+1];
  62   for (Int_t i=0; i<=nBinX; ++i)
  63     binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
  64
  65   Float_t* binLimitsY = new Float_t[nBinY+1];
  66   for (Int_t i=0; i<=nBinY; ++i)
  67     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
  68
  69   Float_t* binLimitsZ = new Float_t[nBinZ+1];
  70   for (Int_t i=0; i<=nBinZ; ++i)
  71     binLimitsZ[i] = Zmin + (Zmax - Zmin) / nBinZ * i;
  72
  73   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
  74
  75   delete[] binLimitsX;
  76   delete[] binLimitsY;
  77   delete[] binLimitsZ;
  78 }
  79
  80 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
  81       Int_t nBinX, Float_t Xmin, Float_t Xmax,
  82       Int_t nBinY, Float_t Ymin, Float_t Ymax,
  83       Int_t nBinZ, const Float_t* zbins)
  84   : AliCorrectionMatrix(name, title)
  85 {
  86   // constructor with variable bin sizes
  87
  88   Float_t* binLimitsX = new Float_t[nBinX+1];
  89   for (Int_t i=0; i<=nBinX; ++i)
  90     binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
  91
  92   Float_t* binLimitsY = new Float_t[nBinY+1];
  93   for (Int_t i=0; i<=nBinY; ++i)
  94     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
  95
  96   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
  97
  98   delete[] binLimitsX;
  99   delete[] binLimitsY;
 100 }
 101
 102 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3F* hBinning)
 103   : AliCorrectionMatrix(name, title)
 104 {
 105   // constructor with variable bin sizes (uses binning of hBinning)
 106
 107   // do not add this hists to the directory
 108   Bool_t oldStatus = TH1::AddDirectoryStatus();
 109   TH1::AddDirectory(kFALSE);
 110
 111   fhMeas = (TH3F*)hBinning->Clone("measured");
 112   fhGene = (TH3F*)hBinning->Clone("generated");
 113   fhCorr = (TH3F*)hBinning->Clone("correction");
 114
 115   fhMeas->SetTitle(Form("%s measured",title)  );
 116   fhGene->SetTitle(Form("%s generated",title ));
 117   fhCorr->SetTitle(Form("%s correction",title));
 118
 119   fhMeas->Reset();
 120   fhGene->Reset();
 121   fhCorr->Reset();
 122
 123   TH1::AddDirectory(oldStatus);
 124
 125   fhMeas->Sumw2();
 126   fhGene->Sumw2();
 127   fhCorr->Sumw2();
 128 }
 129
 130 //____________________________________________________________________
 131 void AliCorrectionMatrix3D::CreateHists(Int_t nBinX, const Float_t* binLimitsX,
 132       Int_t nBinY, const Float_t* binLimitsY,
 133       Int_t nBinZ, const Float_t* binLimitsZ)
 134 {
 135   // create the histograms
 136
 137   // do not add this hists to the directory
 138   Bool_t oldStatus = TH1::AddDirectoryStatus();
 139   TH1::AddDirectory(kFALSE);
 140
 141   fhMeas  = new TH3F("measured",   Form("%s measured",GetTitle()),   nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
 142   fhGene  = new TH3F("generated",  Form("%s generated",GetTitle()),  nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
 143   fhCorr  = new TH3F("correction", Form("%s correction",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
 144
 145   fhMeas->Sumw2();
 146   fhGene->Sumw2();
 147   fhCorr->Sumw2();
 148
 149   TH1::AddDirectory(oldStatus);
 150 }
 151
 152
 153 //____________________________________________________________________
 154 AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
 155 {
 156   //
 157   // destructor
 158   //
 159
 160   // histograms already deleted in base class
 161 }
 162
 163 //____________________________________________________________________
 164 TH3F* AliCorrectionMatrix3D::GetGeneratedHistogram()
 165 {
 166   // return generated histogram casted to correct type
 167   return dynamic_cast<TH3F*> (fhGene);
 168 }
 169
 170 //____________________________________________________________________
 171 TH3F* AliCorrectionMatrix3D::GetMeasuredHistogram()
 172 {
 173   // return measured histogram casted to correct type
 174   return dynamic_cast<TH3F*> (fhMeas);
 175 }
 176
 177 //____________________________________________________________________
 178 TH3F* AliCorrectionMatrix3D::GetCorrectionHistogram()
 179 {
 180   // return correction histogram casted to correct type
 181   return dynamic_cast<TH3F*> (fhCorr);
 182 }
 183
 184 //____________________________________________________________________
 185 AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Char_t* opt, Float_t aMin, Float_t aMax)
 186 {
 187   // returns a 2D projection of this correction
 188
 189   TString option = opt;
 190
 191   if (aMin<aMax) {
 192     if (option.Contains("xy") || option.Contains("yx")) {
 193       Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
 194       Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
 195       fhMeas->GetZaxis()->SetRange(bMin, bMax);
 196     }
 197     else if (option.Contains("xz") || option.Contains("zx")) {
 198       Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
 199       Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
 200       fhMeas->GetYaxis()->SetRange(bMin, bMax);
 201     }
 202     else if (option.Contains("yz") || option.Contains("zy")) {
 203       Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
 204       Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
 205       fhMeas->GetXaxis()->SetRange(bMin, bMax);
 206     }
 207     else {
 208       AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
 209       return 0;
 210     }
 211   }
 212   AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
 213
 214   TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(opt);
 215   TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(opt);
 216
 217   TH2F* corr = (TH2F*)gene->Clone("corr");
 218   corr->Reset();
 219
 220   corr2D->SetGeneratedHistogram(gene);
 221   corr2D->SetMeasuredHistogram(meas);
 222   corr2D->SetCorrectionHistogram(corr);
 223
 224   corr2D->Divide();
 225
 226   // unzoom
 227   fhMeas->GetXaxis()->UnZoom();
 228   fhMeas->GetYaxis()->UnZoom();
 229   fhMeas->GetZaxis()->UnZoom();
 230
 231   fhGene->GetXaxis()->UnZoom();
 232   fhGene->GetYaxis()->UnZoom();
 233   fhGene->GetZaxis()->UnZoom();
 234
 235   return corr2D;
 236 }
 237
 238 //____________________________________________________________________
 239 TH1F* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Char_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
 240 {
 241   // returns a 1D projection of this correction
 242   AliDebug(AliLog::kWarning, Form("WARNING: test"));
 243
 244   AliCorrectionMatrix2D* corr2D;
 245   if (strcmp(opt,"x")==0) {
 246     corr2D = Get2DCorrection("xy",aMin1,aMax1);
 247     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
 248   }
 249   if (strcmp(opt,"y")==0) {
 250     corr2D = Get2DCorrection("xy",aMin1,aMax1);
 251     return corr2D->Get1DCorrectionHistogram("y",aMin2,aMax2);
 252   }
 253   if (strcmp(opt,"z")==0) {
 254     corr2D = Get2DCorrection("yz",aMin1,aMax1);
 255     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
 256   }
 257   AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));
 258
 259   return 0;
 260 }
 261
 262
 263
 264 //____________________________________________________________________
 265 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
 266 {
 267   // add value to measured histogram
 268   GetMeasuredHistogram()->Fill(ax, ay, az);
 269 }
 270
 271 //____________________________________________________________________
 272 void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az)
 273 {
 274   // add value to generated histogram
 275   GetGeneratedHistogram()->Fill(ax, ay, az);
 276 }
 277
 278 //____________________________________________________________________
 279 Float_t AliCorrectionMatrix3D::GetCorrection(Float_t ax, Float_t ay, Float_t az) const
 280 {
 281   // returns a value of the correction map
 282   return fhCorr->GetBinContent(fhCorr->FindBin(ax, ay, az));
 283 }
 284
 285 //____________________________________________________________________
 286 //void AliCorrectionMatrix3D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge, Int_t nBinsZedge)
 287 void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
 288 {
 289   // so what do we do here...
 290 }
 291
 292 //____________________________________________________________________
 293 void AliCorrectionMatrix3D::SaveHistograms()
 294 {
 295   //
 296   // saves the histograms
 297   //
 298
 299   AliCorrectionMatrix::SaveHistograms();
 300
 301   if (GetGeneratedHistogram() && GetMeasuredHistogram())
 302   {
 303     gDirectory->cd(GetName());
 304
 305     AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
 306
 307     gDirectory->cd("..");
 308   }
 309 }
 310
 311 //____________________________________________________________________
 312 Int_t AliCorrectionMatrix3D::CheckEmptyBins(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, Bool_t quiet)
 313 {
 314   //
 315   // counts the number of empty Bins in a given region
 316   //
 317
 318   TH3F* hist = GetGeneratedHistogram();
 319   if (!hist)
 320     return -1;
 321
 322   Int_t emptyBins = 0;
 323   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
 324     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
 325       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
 326         if (hist->GetBinContent(x, y, z) == 0)
 327         {
 328           if (!quiet)
 329             printf("Empty bin in %s at vtx = %f, eta = %f, pt = %f\n", GetName(), hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
 330           ++emptyBins;
 331         }
 332
 333   return emptyBins;
 334 }
 335
 336 //____________________________________________________________________
 337 TH1F* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
 338 {
 339   //
 340   // makes a 1d plots of the relative bin errors
 341   //
 342
 343   TH3F* hist = GetCorrectionHistogram();
 344   if (!hist)
 345     return 0;
 346
 347   TH1F* target = new TH1F("relerrors", "relerrors", 100, 0, 10);
 348
 349   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
 350     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
 351       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
 352         if (hist->GetBinContent(x, y, z) != 0)
 353           target->Fill(100 * hist->GetBinError(x, y, z) / hist->GetBinContent(x, y, z));
 354
 355   return target;
 356 }
 357