]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFECorrectSpectrumBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFECorrectSpectrumBase.cxx
CommitLineData
959ea9d8 1
2/**************************************************************************
3* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4* *
5* Author: The ALICE Off-line Project. *
6* Contributors are mentioned in the code where appropriate. *
7* *
8* Permission to use, copy, modify and distribute this software and its *
9* documentation strictly for non-commercial purposes is hereby granted *
10* without fee, provided that the above copyright notice appears in all *
11* copies and that both the copyright notice and this permission notice *
12* appear in the supporting documentation. The authors make no claims *
13* about the suitability of this software for any purpose. It is *
14* provided "as is" without express or implied warranty. *
15**************************************************************************/
16//
17// Class for spectrum correction
18// Subtraction of hadronic background, Unfolding of the data and
19// Renormalization done here
20// The following containers have to be set:
21// - Correction framework container for real data
22// - Correction framework container for MC (Efficiency Map)
23// - Correction framework container for background coming from data
24// - Correction framework container for background coming from MC
25//
26// Author:
27// Raphaelle Bailhache <R.Bailhache@gsi.de>
28// Markus Fasel <M.Fasel@gsi.de>
29//
30
31#include <TArrayD.h>
32#include <TH1.h>
33#include <TList.h>
34#include <TObjArray.h>
35#include <TROOT.h>
36#include <TCanvas.h>
37#include <TLegend.h>
38#include <TStyle.h>
39#include <TMath.h>
40#include <TAxis.h>
41#include <TGraphErrors.h>
42#include <TFile.h>
43#include <TPad.h>
44#include <TH2D.h>
45#include <TF1.h>
46
47#include "AliPID.h"
48#include "AliCFContainer.h"
49#include "AliCFDataGrid.h"
50#include "AliCFEffGrid.h"
51#include "AliCFGridSparse.h"
52#include "AliCFUnfolding.h"
53#include "AliLog.h"
54
55#include "AliHFECorrectSpectrumBase.h"
56#include "AliHFEcuts.h"
57#include "AliHFEcontainer.h"
58#include "AliHFEtools.h"
59
60ClassImp(AliHFECorrectSpectrumBase)
61
62//____________________________________________________________
63AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const char *name):
64 TNamed(name, ""),
65 fCFContainers(new TObjArray(kDataContainerV0+1)),
66 fCorrelation(NULL),
67 fEfficiencyFunction(NULL),
68 fEtaSelected(kFALSE),
69 fSetSmoothing(kFALSE),
70 fNbDimensions(1),
71 fNEvents(0),
72 fStepMC(-1),
73 fStepTrue(-1),
74 fStepData(-1),
75 fStepBeforeCutsV0(-1),
76 fStepAfterCutsV0(-1),
77 fStepGuessedUnfolding(-1),
78 fNumberOfIterations(10),
79 fChargeChoosen(kAllCharge),
80 fTestCentralityLow(-1),
81 fTestCentralityHigh(-1)
82{
83 //
84 // Default constructor
85 //
86
87 memset(fEtaRange, 0, sizeof(Double_t) * 2);
88 memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
89
90}
91//____________________________________________________________
92AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref):
93 TNamed(ref),
94 fCFContainers(NULL),
95 fCorrelation(NULL),
96 fEfficiencyFunction(NULL),
97 fEtaSelected(ref.fEtaSelected),
98 fSetSmoothing(ref.fSetSmoothing),
99 fNbDimensions(ref.fNbDimensions),
100 fNEvents(ref.fNEvents),
101 fStepMC(ref.fStepMC),
102 fStepTrue(ref.fStepTrue),
103 fStepData(ref.fStepData),
104 fStepBeforeCutsV0(ref.fStepBeforeCutsV0),
105 fStepAfterCutsV0(ref.fStepAfterCutsV0),
106 fStepGuessedUnfolding(ref.fStepGuessedUnfolding),
107 fNumberOfIterations(ref.fNumberOfIterations),
108 fChargeChoosen(ref.fChargeChoosen),
109 fTestCentralityLow(ref.fTestCentralityLow),
110 fTestCentralityHigh(ref.fTestCentralityHigh)
111{
112 //
113 // Copy constructor
114 //
115 ref.Copy(*this);
116
117}
118//____________________________________________________________
119AliHFECorrectSpectrumBase &AliHFECorrectSpectrumBase::operator=(const AliHFECorrectSpectrumBase &ref){
120 //
121 // Assignment operator
122 //
123 if(this == &ref)
124 ref.Copy(*this);
125 return *this;
126}
127//____________________________________________________________
128void AliHFECorrectSpectrumBase::Copy(TObject &o) const {
129 //
130 // Copy into object o
131 //
132 AliHFECorrectSpectrumBase &target = dynamic_cast<AliHFECorrectSpectrumBase &>(o);
133 target.fCFContainers = fCFContainers;
134 target.fCorrelation = fCorrelation;
135 target.fEfficiencyFunction = fEfficiencyFunction;
136 target.fEtaSelected = fEtaSelected;
137 target.fSetSmoothing = fSetSmoothing;
138 target.fNbDimensions = fNbDimensions;
139 target.fNEvents = fNEvents;
140 target.fStepMC = fStepMC;
141 target.fStepTrue = fStepTrue;
142 target.fStepData = fStepData;
143 target.fStepBeforeCutsV0 = fStepBeforeCutsV0;
144 target.fStepAfterCutsV0 = fStepAfterCutsV0;
145 target.fStepGuessedUnfolding = fStepGuessedUnfolding;
146 target.fNumberOfIterations = fNumberOfIterations;
147 target.fChargeChoosen = fChargeChoosen;
148 target.fTestCentralityLow = fTestCentralityLow;
149 target.fTestCentralityHigh = fTestCentralityHigh;
150 target.fEtaRange[0] = fEtaRange[0];
151 target.fEtaRange[1] = fEtaRange[1];
152 target.fEtaRangeNorm[0] = fEtaRangeNorm[0];
153 target.fEtaRangeNorm[1] = fEtaRangeNorm[1];
154
155}
156
157//____________________________________________________________
158AliHFECorrectSpectrumBase::~AliHFECorrectSpectrumBase(){
159 //
160 // Destructor
161 //
162 if(fCFContainers) delete fCFContainers;
163
164}
165//__________________________________________________________________________________
166TGraphErrors *AliHFECorrectSpectrumBase::Normalize(THnSparse * const spectrum) const {
167 //
168 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
169 // Give the final pt spectrum to be compared
170 //
171
172 if(fNEvents > 0) {
173
174 TH1D* projection = spectrum->Projection(0);
175 CorrectFromTheWidth(projection);
176 TGraphErrors *graphError = NormalizeTH1(projection);
177 return graphError;
178
179 }
180
181 return 0x0;
182
183
184}
185//__________________________________________________________________________________
186TGraphErrors *AliHFECorrectSpectrumBase::Normalize(AliCFDataGrid * const spectrum) const {
187 //
188 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
189 // Give the final pt spectrum to be compared
190 //
191 if(fNEvents > 0) {
192
193 TH1D* projection = (TH1D *) spectrum->Project(0);
194 CorrectFromTheWidth(projection);
195 TGraphErrors *graphError = NormalizeTH1(projection);
196
197 return graphError;
198
199 }
200
201 return 0x0;
202
203
204}
205//__________________________________________________________________________________
206TGraphErrors *AliHFECorrectSpectrumBase::NormalizeTH1(TH1 *input) const {
207 //
208 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
209 // Give the final pt spectrum to be compared
210 //
211 Double_t chargecoefficient = 0.5;
212 if(fChargeChoosen != 0) chargecoefficient = 1.0;
213
214 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
215 printf("Normalizing Eta Range %f\n", etarange);
216 if(fNEvents > 0) {
217
218 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
219 Double_t p = 0, dp = 0; Int_t point = 1;
220 Double_t n = 0, dN = 0;
221 Double_t nCorr = 0, dNcorr = 0;
222 Double_t errdN = 0, errdp = 0;
223 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
224 point = ibin - input->GetXaxis()->GetFirst();
225 p = input->GetXaxis()->GetBinCenter(ibin);
226 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
227 n = input->GetBinContent(ibin);
228 dN = input->GetBinError(ibin);
229 // New point
230 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * 1./(2. * TMath::Pi() * p) * n;
231 errdN = 1./(2. * TMath::Pi() * p);
232 errdp = 1./(2. * TMath::Pi() * p*p) * n;
233 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
234
235 spectrumNormalized->SetPoint(point, p, nCorr);
236 spectrumNormalized->SetPointError(point, dp, dNcorr);
237 }
238 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
239 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
240 spectrumNormalized->SetMarkerStyle(22);
241 spectrumNormalized->SetMarkerColor(kBlue);
242 spectrumNormalized->SetLineColor(kBlue);
243
244 return spectrumNormalized;
245
246 }
247
248 return 0x0;
249
250
251}
252//____________________________________________________________
253void AliHFECorrectSpectrumBase::SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type){
254 //
255 // Set the container for a given step to the
256 //
257 if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
258 fCFContainers->AddAt(cont, type);
259}
260
261//____________________________________________________________
262AliCFContainer *AliHFECorrectSpectrumBase::GetContainer(AliHFECorrectSpectrumBase::CFContainer_t type){
263 //
264 // Get Correction Framework Container for given type
265 //
266 if(!fCFContainers) return NULL;
267 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
268}
269//____________________________________________________________
270AliCFContainer *AliHFECorrectSpectrumBase::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
271 //
272 // Slice bin for a given source of electron
273 // nDim is the number of dimension the corrections are done
274 // dimensions are the definition of the dimensions
275 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
276 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
277 // centrality (-1 means we do not cut on centrality)
278 //
279
280 Double_t *varMin = new Double_t[container->GetNVar()],
281 *varMax = new Double_t[container->GetNVar()];
282
283 Double_t *binLimits;
284 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
285
286 binLimits = new Double_t[container->GetNBins(ivar)+1];
287 container->GetBinLimits(ivar,binLimits);
288 varMin[ivar] = binLimits[0];
289 varMax[ivar] = binLimits[container->GetNBins(ivar)];
290 // source
291 if(ivar == 4){
292 if((source>= 0) && (source<container->GetNBins(ivar))) {
293 varMin[ivar] = binLimits[source];
294 varMax[ivar] = binLimits[source];
295 }
296 }
297 // charge
298 if(ivar == 3) {
299 if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
300 }
301 // eta
302 if(ivar == 1){
303 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
304 AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
305 if(fEtaSelected){
306 varMin[ivar] = fEtaRange[0];
307 varMax[ivar] = fEtaRange[1];
308 }
309 }
310 if(fEtaSelected){
311 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
312 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
313 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
314 }
315 // centrality
316 if(ivar == 5){
317 if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
318 varMin[ivar] = binLimits[centralitylow];
319 varMax[ivar] = binLimits[centralityhigh];
320
321 TAxis *axistest = container->GetAxis(5,0);
322 AliDebug(1, Form("Number of bin in centrality direction %d\n",axistest->GetNbins()));
323 AliDebug(1, Form("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]));
324 Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
325 Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
326 AliDebug(1, Form("Low centrality %f and high centrality %f\n",lowcentrality,highcentrality));
327
328 }
329 }
330
331 delete[] binLimits;
332
333 }
334
335 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
336 delete[] varMin; delete[] varMax;
337
338 return k;
339
340}
341
342//_________________________________________________________________________
343THnSparseF *AliHFECorrectSpectrumBase::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions,Chargetype_t charge,Int_t centralitylow, Int_t centralityhigh) const {
344 //
345 // Slice correlation
346 //
347
348 Int_t ndimensions = correlationmatrix->GetNdimensions();
349 //printf("Number of dimension %d correlation map\n",ndimensions);
350 if(ndimensions < (2*nDim)) {
351 AliError("Problem in the dimensions");
352 return NULL;
353 }
354
355 // Cut in centrality is centrality > -1
356 if((5+((Int_t)(ndimensions/2.))) < ndimensions) {
357 if((centralitylow >=0) && (centralityhigh >=0)) {
358
359 TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
360 TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
361
362 Int_t bins0 = axiscentrality0->GetNbins();
363 Int_t bins1 = axiscentrality1->GetNbins();
364
365 AliDebug(1, Form("Number of centrality bins: %d and %d\n",bins0,bins1));
366 if(bins0 != bins1) {
367 AliError("Problem in the dimensions");
368 return NULL;
369 }
370
371 if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
372 axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
373 axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
374
375 Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
376 Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
377 Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
378 Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
379 AliDebug(1,Form("0 Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0));
380 AliDebug(1,Form("1 Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1));
381
382 }
383 }
384 }
385
386 // Cut in eta > -1
387 if(fEtaSelected){
388 if((1+((Int_t)(ndimensions/2.))) < ndimensions) {
389
390 TAxis *axiseta0 = correlationmatrix->GetAxis(1);
391 TAxis *axiseta1 = correlationmatrix->GetAxis(1+((Int_t)(ndimensions/2.)));
392
393 Int_t bins0 = axiseta0->GetNbins();
394 Int_t bins1 = axiseta1->GetNbins();
395
396 AliDebug(1, Form("Number of eta bins: %d and %d\n",bins0,bins1));
397 if(bins0 != bins1) {
398 AliError("Problem in the dimensions");
399 return NULL;
400 }
401
402 axiseta0->SetRangeUser(fEtaRange[0],fEtaRange[1]);
403 axiseta1->SetRangeUser(fEtaRange[0],fEtaRange[1]);
404
405 Double_t loweta0 = axiseta0->GetBinLowEdge(axiseta0->FindBin(fEtaRange[0]));
406 Double_t higheta0 = axiseta0->GetBinUpEdge(axiseta0->FindBin(fEtaRange[1]));
407 Double_t loweta1 = axiseta1->GetBinLowEdge(axiseta1->FindBin(fEtaRange[0]));
408 Double_t higheta1 = axiseta1->GetBinUpEdge(axiseta1->FindBin(fEtaRange[1]));
409 AliInfo(Form("0 Low eta %f and high eta %f\n",loweta0,higheta0));
410 AliInfo(Form("1 Low eta %f and high eta %f\n",loweta1,higheta1));
411
412 }
413 }
414
415 // Cut in charge
416 if(charge != kAllCharge) {
417 if((3+((Int_t)(ndimensions/2.))) < ndimensions) {
418
419 TAxis *axischarge0 = correlationmatrix->GetAxis(3);
420 TAxis *axischarge1 = correlationmatrix->GetAxis(3+((Int_t)(ndimensions/2.)));
421
422 Int_t bins0 = axischarge0->GetNbins();
423 Int_t bins1 = axischarge1->GetNbins();
424
425 AliDebug(1, Form("Number of charge bins: %d and %d\n",bins0,bins1));
426 if(bins0 != bins1) {
427 AliError("Problem in the dimensions");
428 return NULL;
429 }
430
431 axischarge0->SetRangeUser(charge,charge);
432 axischarge1->SetRangeUser(charge,charge);
433
434 Double_t lowcharge0 = axischarge0->GetBinLowEdge(axischarge0->FindBin(charge));
435 Double_t highcharge0 = axischarge0->GetBinUpEdge(axischarge0->FindBin(charge));
436 Double_t lowcharge1 = axischarge1->GetBinLowEdge(axischarge1->FindBin(charge));
437 Double_t highcharge1 = axischarge1->GetBinUpEdge(axischarge1->FindBin(charge));
438 AliInfo(Form("0 Low charge %f and high charge %f\n",lowcharge0,highcharge0));
439 AliInfo(Form("1 Low charge %f and high charge %f\n",lowcharge1,highcharge1));
440
441 }
442 }
443
444
445 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
446
447 Int_t *dim = new Int_t[nDim*2];
448 for(Int_t iter=0; iter < nDim; iter++){
449 dim[iter] = dimensions[iter];
450 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
451 }
452
453 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
454
455 delete[] dim;
456 return k;
457
458}
459//___________________________________________________________________________
460void AliHFECorrectSpectrumBase::CorrectFromTheWidth(TH1D *h1) const {
461 //
462 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
463 //
464
465 TAxis *axis = h1->GetXaxis();
466 Int_t nbinX = h1->GetNbinsX();
467
468 for(Int_t i = 1; i <= nbinX; i++) {
469
470 Double_t width = axis->GetBinWidth(i);
471 Double_t content = h1->GetBinContent(i);
472 Double_t error = h1->GetBinError(i);
473 h1->SetBinContent(i,content/width);
474 h1->SetBinError(i,error/width);
475 }
476
477}
478
479//___________________________________________________________________________
480void AliHFECorrectSpectrumBase::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
481 //
482 // Correct statistical error
483 //
484
485 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
486 Int_t nbinX = h1->GetNbinsX();
487 Int_t bins[1];
488 for(Long_t i = 1; i <= nbinX; i++) {
489 bins[0] = i;
490 Float_t content = h1->GetBinContent(i);
491 if(content>0){
492 Float_t error = TMath::Sqrt(content);
493 backgroundGrid->SetElementError(bins, error);
494 }
495 }
496}
497//_________________________________________________________________________
498TObject* AliHFECorrectSpectrumBase::GetSpectrum(const AliCFContainer * const c, Int_t step) {
499 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
500 return data;
501}
502//_________________________________________________________________________
503TObject* AliHFECorrectSpectrumBase::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
504 //
505 // Create efficiency grid and calculate efficiency
506 // of step to step0
507 //
508 TString name("eff");
509 name += step;
510 name+= step0;
511 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
512 eff->CalculateEfficiency(step,step0);
513 return eff;
514}