PWG0/AliCorrectionMatrix3D.cxx

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