New set of default cuts suitable for pp events.
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix3D.cxx
1 /* $Id$ */
2
3 // ------------------------------------------------------
4 //
5 // Class to handle 3d-corrections.
6 //
7 // ------------------------------------------------------
8 //
9
10 #include <TH3F.h>
11 #include <TH1F.h>
12
13 #include <AliLog.h>
14
15 #include "AliCorrectionMatrix3D.h"
16 #include "AliPWG0Helper.h"
17
18 //____________________________________________________________________
19 ClassImp(AliCorrectionMatrix3D)
20
21 //____________________________________________________________________
22 AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
23   AliCorrectionMatrix()
24 {
25   // default constructor
26 }
27
28 //____________________________________________________________________
29 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
30   : AliCorrectionMatrix(c)
31 {
32   // copy constructor
33   ((AliCorrectionMatrix3D &)c).Copy(*this);
34 }
35
36 //____________________________________________________________________
37 AliCorrectionMatrix3D &AliCorrectionMatrix3D::operator=(const AliCorrectionMatrix3D &c)
38 {
39   // assigment operator
40
41   if (this != &c)
42     ((AliCorrectionMatrix3D &) c).Copy(*this);
43
44   return *this;
45 }
46
47 //____________________________________________________________________
48 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
49               Int_t nBinX, Float_t Xmin, Float_t Xmax,
50               Int_t nBinY, Float_t Ymin, Float_t Ymax,
51               Int_t nBinZ, Float_t Zmin, Float_t Zmax)
52   : AliCorrectionMatrix(name, title)
53 {
54   //
55   // constructor
56   //
57
58   // do not add this hists to the directory
59   Bool_t oldStatus = TH1::AddDirectoryStatus();
60   TH1::AddDirectory(kFALSE);
61
62   fhMeas  = new TH3F("measured",   Form("%s measured",title),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
63   fhGene  = new TH3F("generated",  Form("%s generated",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
64   fhCorr  = new TH3F("correction", Form("%s correction",title), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax, nBinZ, Zmin, Zmax);
65
66   TH1::AddDirectory(oldStatus);
67
68   fhMeas->Sumw2();
69   fhGene->Sumw2();
70   fhCorr->Sumw2();
71 }
72
73 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
74       Int_t nBinX, Float_t Xmin, Float_t Xmax,
75       Int_t nBinY, Float_t Ymin, Float_t Ymax,
76       Int_t nBinZ, const Float_t* zbins)
77   : AliCorrectionMatrix(name, title)
78 {
79   // constructor with variable bin sizes
80
81   Float_t* binLimitsX = new Float_t[nBinX+1];
82   for (Int_t i=0; i<=nBinX; ++i)
83     binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
84
85   Float_t* binLimitsY = new Float_t[nBinY+1];
86   for (Int_t i=0; i<=nBinY; ++i)
87     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
88
89   // do not add this hists to the directory
90   Bool_t oldStatus = TH1::AddDirectoryStatus();
91   TH1::AddDirectory(kFALSE);
92
93   fhMeas  = new TH3F("measured",   Form("%s measured",title),   nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
94   fhGene  = new TH3F("generated",  Form("%s generated",title),  nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
95   fhCorr  = new TH3F("correction", Form("%s correction",title), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
96
97   TH1::AddDirectory(oldStatus);
98
99   delete[] binLimitsX;
100   delete[] binLimitsY;
101
102   fhMeas->Sumw2();
103   fhGene->Sumw2();
104   fhCorr->Sumw2();
105 }
106
107 //____________________________________________________________________
108 AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
109 {
110   //
111   // destructor
112   //
113
114   // histograms already deleted in base class
115 }
116
117 //____________________________________________________________________
118 TH3F* AliCorrectionMatrix3D::GetGeneratedHistogram()
119 {
120   // return generated histogram casted to correct type
121   return dynamic_cast<TH3F*> (fhGene);
122 }
123
124 //____________________________________________________________________
125 TH3F* AliCorrectionMatrix3D::GetMeasuredHistogram()
126 {
127   // return measured histogram casted to correct type
128   return dynamic_cast<TH3F*> (fhMeas);
129 }
130
131 //____________________________________________________________________
132 TH3F* AliCorrectionMatrix3D::GetCorrectionHistogram()
133 {
134   // return correction histogram casted to correct type
135   return dynamic_cast<TH3F*> (fhCorr);
136 }
137
138 //____________________________________________________________________
139 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
140 {
141   // add value to measured histogram
142   GetMeasuredHistogram()->Fill(ax, ay, az);
143 }
144
145 //____________________________________________________________________
146 void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az)
147 {
148   // add value to generated histogram
149   GetGeneratedHistogram()->Fill(ax, ay, az);
150 }
151
152 //____________________________________________________________________
153 Float_t AliCorrectionMatrix3D::GetCorrection(Float_t ax, Float_t ay, Float_t az) const
154 {
155   // returns a value of the correction map
156   return fhCorr->GetBinContent(fhCorr->FindBin(ax, ay, az));
157 }
158
159 //____________________________________________________________________
160 //void AliCorrectionMatrix3D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge, Int_t nBinsZedge)
161 void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
162 {
163   // so what do we do here...
164 }
165
166 //____________________________________________________________________
167 void AliCorrectionMatrix3D::SaveHistograms()
168 {
169   //
170   // saves the histograms
171   //
172
173   AliCorrectionMatrix::SaveHistograms();
174
175   if (GetGeneratedHistogram() && GetMeasuredHistogram())
176   {
177     gDirectory->cd(GetName());
178     
179     AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
180     
181     gDirectory->cd("..");
182   }
183 }
184
185 //____________________________________________________________________
186 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)
187 {
188   //
189   // counts the number of empty Bins in a given region
190   //
191
192   TH3F* hist = GetGeneratedHistogram();
193   if (!hist)
194     return -1;
195
196   Int_t emptyBins = 0;
197   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
198     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
199       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
200         if (hist->GetBinContent(x, y, z) == 0)
201         {
202           if (!quiet)
203             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));
204           ++emptyBins;
205         }
206
207   return emptyBins;
208 }
209
210 //____________________________________________________________________
211 TH1F* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
212 {
213   //
214   // makes a 1d plots of the relative bin errors
215   //
216
217   TH3F* hist = GetCorrectionHistogram();
218   if (!hist)
219     return 0;
220     
221   TH1F* target = new TH1F("relerrors", "relerrors", 100, 0, 10);
222
223   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
224     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
225       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
226         if (hist->GetBinContent(x, y, z) != 0)
227           target->Fill(100 * hist->GetBinError(x, y, z) / hist->GetBinContent(x, y, z));
228
229   return target;
230 }
231