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