]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/AliCorrectionMatrix3D.cxx
adding script to retrieve all datapoints for testing purposes
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix3D.cxx
index a714603956c4a4352d96a41c4b03b0ede7b9bde3..b953f3b89909cece7c5b4e8edee3de424c8a32b8 100644 (file)
@@ -7,12 +7,16 @@
 // ------------------------------------------------------
 //
 
-#include <TH3F.h>
+#include <TDirectory.h>
 #include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TString.h>
 
 #include <AliLog.h>
 
 #include "AliCorrectionMatrix3D.h"
+#include "AliCorrectionMatrix2D.h"
 #include "AliPWG0Helper.h"
 
 //____________________________________________________________________
@@ -55,13 +59,23 @@ AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* t
   // constructor
   //
 
-  fhMeas  = new TH3F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
-  fhGene  = new TH3F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
-  fhCorr  = new TH3F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
+  Float_t* binLimitsX = new Float_t[nBinX+1];
+  for (Int_t i=0; i<=nBinX; ++i)
+    binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
 
-  fhMeas->Sumw2();
-  fhGene->Sumw2();
-  fhCorr->Sumw2();
+  Float_t* binLimitsY = new Float_t[nBinY+1];
+  for (Int_t i=0; i<=nBinY; ++i)
+    binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
+
+  Float_t* binLimitsZ = new Float_t[nBinZ+1];
+  for (Int_t i=0; i<=nBinZ; ++i)
+    binLimitsZ[i] = Zmin + (Zmax - Zmin) / nBinZ * i;
+
+  CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
+
+  delete[] binLimitsX;
+  delete[] binLimitsY;
+  delete[] binLimitsZ;
 }
 
 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
@@ -80,18 +94,63 @@ AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* t
   for (Int_t i=0; i<=nBinY; ++i)
     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
 
-  fhMeas  = new TH3F(Form("meas_%s",name), Form("meas_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
-  fhGene  = new TH3F(Form("gene_%s",name), Form("gene_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
-  fhCorr  = new TH3F(Form("corr_%s",name), Form("corr_%s",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
+  CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
 
   delete[] binLimitsX;
   delete[] binLimitsY;
+}
+
+AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3F* hBinning)
+  : AliCorrectionMatrix(name, title)
+{
+  // constructor with variable bin sizes (uses binning of hBinning)
+
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  fhMeas = (TH3F*)hBinning->Clone("measured");
+  fhGene = (TH3F*)hBinning->Clone("generated");
+  fhCorr = (TH3F*)hBinning->Clone("correction");
+
+  fhMeas->SetTitle(Form("%s measured", GetTitle()));
+  fhGene->SetTitle(Form("%s generated", GetTitle()));
+  fhCorr->SetTitle(Form("%s correction", GetTitle()));
+
+  fhMeas->Reset();
+  fhGene->Reset();
+  fhCorr->Reset();
+
+  TH1::AddDirectory(oldStatus);
 
   fhMeas->Sumw2();
   fhGene->Sumw2();
   fhCorr->Sumw2();
 }
 
+//____________________________________________________________________
+void AliCorrectionMatrix3D::CreateHists(Int_t nBinX, const Float_t* binLimitsX,
+      Int_t nBinY, const Float_t* binLimitsY,
+      Int_t nBinZ, const Float_t* binLimitsZ)
+{
+  // create the histograms
+
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  fhMeas  = new TH3F("measured",   Form("%s measured",GetTitle()),   nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
+  fhGene  = new TH3F("generated",  Form("%s generated",GetTitle()),  nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
+  fhCorr  = new TH3F("correction", Form("%s correction",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
+
+  fhMeas->Sumw2();
+  fhGene->Sumw2();
+  fhCorr->Sumw2();
+
+  TH1::AddDirectory(oldStatus);
+}
+
+
 //____________________________________________________________________
 AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
 {
@@ -123,6 +182,97 @@ TH3F* AliCorrectionMatrix3D::GetCorrectionHistogram()
   return dynamic_cast<TH3F*> (fhCorr);
 }
 
+//____________________________________________________________________
+AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Char_t* opt, Float_t aMin, Float_t aMax)
+{
+  // returns a 2D projection of this correction
+
+  TString option = opt;
+
+  // unzoom
+  fhMeas->GetXaxis()->UnZoom();
+  fhMeas->GetYaxis()->UnZoom();
+  fhMeas->GetZaxis()->UnZoom();
+
+  fhGene->GetXaxis()->UnZoom();
+  fhGene->GetYaxis()->UnZoom();
+  fhGene->GetZaxis()->UnZoom();
+
+  if (aMin<aMax) {
+    if (option.Contains("xy") || option.Contains("yx")) {
+      Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
+      Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
+      fhGene->GetZaxis()->SetRange(bMin, bMax);
+      fhMeas->GetZaxis()->SetRange(bMin, bMax);
+    }
+    else if (option.Contains("xz") || option.Contains("zx")) {
+      Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
+      Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
+      fhGene->GetYaxis()->SetRange(bMin, bMax);
+      fhMeas->GetYaxis()->SetRange(bMin, bMax);
+    }
+    else if (option.Contains("yz") || option.Contains("zy")) {
+      Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
+      Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
+      fhGene->GetXaxis()->SetRange(bMin, bMax);
+      fhMeas->GetXaxis()->SetRange(bMin, bMax);
+    }
+    else {
+      AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
+      return 0;
+    }
+  }
+
+  AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
+
+  TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(option)->Clone(Form("%s_meas", corr2D->GetName()));
+  TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(option)->Clone(Form("%s_gene", corr2D->GetName()));
+
+  TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
+  corr->Reset();
+
+  corr2D->SetGeneratedHistogram(gene);
+  corr2D->SetMeasuredHistogram(meas);
+  corr2D->SetCorrectionHistogram(corr);
+
+  corr2D->Divide();
+
+  // unzoom
+  fhMeas->GetXaxis()->UnZoom();
+  fhMeas->GetYaxis()->UnZoom();
+  fhMeas->GetZaxis()->UnZoom();
+
+  fhGene->GetXaxis()->UnZoom();
+  fhGene->GetYaxis()->UnZoom();
+  fhGene->GetZaxis()->UnZoom();
+
+  return corr2D;
+}
+
+//____________________________________________________________________
+TH1F* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Char_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
+{
+  // returns a 1D projection of this correction
+  AliDebug(AliLog::kWarning, Form("WARNING: test"));
+  
+  AliCorrectionMatrix2D* corr2D;
+  if (strcmp(opt,"x")==0) {
+    corr2D = Get2DCorrection("yx",aMin1,aMax1);
+    return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
+  }
+  if (strcmp(opt,"y")==0) {
+    corr2D = Get2DCorrection("xy",aMin1,aMax1);
+    return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
+  }  
+  if (strcmp(opt,"z")==0) {
+    corr2D = Get2DCorrection("yz",aMin1,aMax1);
+    return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
+  }  
+  AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));  
+  
+  return 0;
+}
+
 //____________________________________________________________________
 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
 {
@@ -161,7 +311,13 @@ void AliCorrectionMatrix3D::SaveHistograms()
   AliCorrectionMatrix::SaveHistograms();
 
   if (GetGeneratedHistogram() && GetMeasuredHistogram())
-    AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram());
+  {
+    gDirectory->cd(GetName());
+    
+    AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
+    
+    gDirectory->cd("..");
+  }
 }
 
 //____________________________________________________________________