2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
27 // Raphaelle Bailhache <R.Bailhache@gsi.de>
28 // Markus Fasel <M.Fasel@gsi.de>
34 #include <TObjArray.h>
41 #include <TGraphErrors.h>
48 #include "AliCFContainer.h"
49 #include "AliCFDataGrid.h"
50 #include "AliCFEffGrid.h"
51 #include "AliCFGridSparse.h"
52 #include "AliCFUnfolding.h"
55 #include "AliHFECorrectSpectrumBase.h"
56 #include "AliHFEcuts.h"
57 #include "AliHFEcontainer.h"
58 #include "AliHFEtools.h"
60 ClassImp(AliHFECorrectSpectrumBase)
62 //____________________________________________________________
63 AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const char *name):
65 fCFContainers(new TObjArray(kNbCFContainers)),
67 fEfficiencyFunction(NULL),
69 fSetSmoothing(kFALSE),
75 fStepBeforeCutsV0(-1),
77 fStepGuessedUnfolding(-1),
78 fNumberOfIterations(10),
79 fChargeChoosen(kAllCharge),
80 fTestCentralityLow(-1),
81 fTestCentralityHigh(-1)
84 // Default constructor
87 memset(fEtaRange, 0, sizeof(Double_t) * 2);
88 memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
91 //____________________________________________________________
92 AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref):
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)
118 //____________________________________________________________
119 AliHFECorrectSpectrumBase &AliHFECorrectSpectrumBase::operator=(const AliHFECorrectSpectrumBase &ref){
121 // Assignment operator
127 //____________________________________________________________
128 void AliHFECorrectSpectrumBase::Copy(TObject &o) const {
130 // Copy into object o
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];
157 //____________________________________________________________
158 AliHFECorrectSpectrumBase::~AliHFECorrectSpectrumBase(){
162 if(fCFContainers) delete fCFContainers;
165 //__________________________________________________________________________________
166 TGraphErrors *AliHFECorrectSpectrumBase::Normalize(THnSparse * const spectrum) const {
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
174 TH1D* projection = spectrum->Projection(0);
175 CorrectFromTheWidth(projection);
176 TGraphErrors *graphError = NormalizeTH1(projection);
185 //__________________________________________________________________________________
186 TGraphErrors *AliHFECorrectSpectrumBase::Normalize(AliCFDataGrid * const spectrum) const {
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
193 TH1D* projection = (TH1D *) spectrum->Project(0);
194 CorrectFromTheWidth(projection);
195 TGraphErrors *graphError = NormalizeTH1(projection);
205 //__________________________________________________________________________________
206 TGraphErrors *AliHFECorrectSpectrumBase::NormalizeTH1(TH1 *input) const {
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
211 Double_t chargecoefficient = 0.5;
212 if(fChargeChoosen != 0) chargecoefficient = 1.0;
214 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
215 printf("Normalizing Eta Range %f\n", etarange);
216 printf("Number of events in Normalisation: %d\n", fNEvents);
217 AliDebug(3, Form("charge coefficient: %f\n", chargecoefficient));
220 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
221 Double_t p = 0, dp = 0; Int_t point = 1;
222 Double_t n = 0, dN = 0;
223 Double_t nCorr = 0, dNcorr = 0;
224 //Double_t errdN = 0, errdp = 0;
226 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
227 point = ibin - input->GetXaxis()->GetFirst();
228 p = input->GetXaxis()->GetBinCenter(ibin);
229 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
230 n = input->GetBinContent(ibin);
231 AliDebug(6, Form("p: %f, n: %e\n", p, n));
232 dN = input->GetBinError(ibin);
234 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * 1./(2. * TMath::Pi() * p) * n;
235 errdN = 1./(2. * TMath::Pi() * p);
236 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
237 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
238 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN);
240 spectrumNormalized->SetPoint(point, p, nCorr);
241 spectrumNormalized->SetPointError(point, dp, dNcorr);
243 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
244 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
245 spectrumNormalized->SetMarkerStyle(22);
246 spectrumNormalized->SetMarkerColor(kBlue);
247 spectrumNormalized->SetLineColor(kBlue);
249 return spectrumNormalized;
257 //____________________________________________________________
258 void AliHFECorrectSpectrumBase::SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type){
260 // Set the container for a given step to the
262 if(!fCFContainers) fCFContainers = new TObjArray(kNbCFContainers);
263 fCFContainers->AddAt(cont, type);
266 //____________________________________________________________
267 AliCFContainer *AliHFECorrectSpectrumBase::GetContainer(AliHFECorrectSpectrumBase::CFContainer_t type){
269 // Get Correction Framework Container for given type
271 if(!fCFContainers) return NULL;
272 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
274 //____________________________________________________________
275 AliCFContainer *AliHFECorrectSpectrumBase::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) {
277 // Slice bin for a given source of electron
278 // nDim is the number of dimension the corrections are done
279 // dimensions are the definition of the dimensions
280 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
281 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
282 // centrality (-1 means we do not cut on centrality)
284 const double kVerySmall = 1e-5;
286 Double_t *varMin = new Double_t[container->GetNVar()],
287 *varMax = new Double_t[container->GetNVar()];
290 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
292 binLimits = new Double_t[container->GetNBins(ivar)+1];
293 container->GetBinLimits(ivar,binLimits);
294 varMin[ivar] = binLimits[0];
295 varMax[ivar] = binLimits[container->GetNBins(ivar)];
298 if((source>= 0) && (source<container->GetNBins(ivar))) {
299 varMin[ivar] = binLimits[source];
300 varMax[ivar] = binLimits[source];
305 if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
309 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
310 AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
312 varMin[ivar] = fEtaRange[0];
313 varMax[ivar] = fEtaRange[1];
317 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
318 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
319 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[1]));
323 if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar)) && doCentralityProjection) {
324 varMin[ivar] = binLimits[centralitylow]+ kVerySmall;
325 varMax[ivar] = binLimits[centralityhigh]+ kVerySmall;
327 TAxis *axistest = container->GetAxis(5,0);
328 AliDebug(1, Form("Number of bin in centrality direction %d\n",axistest->GetNbins()));
329 AliDebug(1, Form("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]));
330 Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
331 Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
332 AliDebug(1, Form("Low centrality %f and high centrality %f\n",lowcentrality,highcentrality));
336 // Protect varmax against overflow
337 if(TMath::Abs(varMax[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
338 AliInfo("Protection against overflow bin");
339 varMax[ivar] -= kVerySmall;
341 if(varMax[ivar] > binLimits[container->GetNBins(ivar)]){
342 AliError("Upper limit exceeds allowed range");
343 varMax[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
346 // Protect varmin against overflow
347 if(TMath::Abs(varMin[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
348 AliInfo("Protection against overflow bin");
349 varMin[ivar] -= kVerySmall;
351 if(varMin[ivar] > binLimits[container->GetNBins(ivar)]){
352 AliError("Upper limit exceeds allowed range");
353 varMin[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
355 AliDebug(1, Form("variable %d: Settting limits to %f and %f\n", ivar, varMin[ivar], varMax[ivar]));
360 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
361 delete[] varMin; delete[] varMax;
367 //_________________________________________________________________________
368 THnSparseF *AliHFECorrectSpectrumBase::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions,Chargetype_t charge,Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) const {
373 Int_t ndimensions = correlationmatrix->GetNdimensions();
374 //printf("Number of dimension %d correlation map\n",ndimensions);
375 if(ndimensions < (2*nDim)) {
376 AliError("Problem in the dimensions");
380 // Cut in centrality is centrality > -1
381 if((5+((Int_t)(ndimensions/2.))) < ndimensions) {
382 if((centralitylow >=0) && (centralityhigh >=0)) {
384 TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
385 TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
387 Int_t bins0 = axiscentrality0->GetNbins();
388 Int_t bins1 = axiscentrality1->GetNbins();
390 AliDebug(1, Form("Number of centrality bins: %d and %d\n",bins0,bins1));
392 AliError("Problem in the dimensions");
396 if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0) && doCentralityProjection) {
397 axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
398 axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
400 Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
401 Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
402 Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
403 Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
404 AliDebug(1,Form("0 Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0));
405 AliDebug(1,Form("1 Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1));
413 if((1+((Int_t)(ndimensions/2.))) < ndimensions) {
415 TAxis *axiseta0 = correlationmatrix->GetAxis(1);
416 TAxis *axiseta1 = correlationmatrix->GetAxis(1+((Int_t)(ndimensions/2.)));
418 Int_t bins0 = axiseta0->GetNbins();
419 Int_t bins1 = axiseta1->GetNbins();
421 AliDebug(1, Form("Number of eta bins: %d and %d\n",bins0,bins1));
423 AliError("Problem in the dimensions");
427 axiseta0->SetRangeUser(fEtaRange[0],fEtaRange[1]);
428 axiseta1->SetRangeUser(fEtaRange[0],fEtaRange[1]);
430 Double_t loweta0 = axiseta0->GetBinLowEdge(axiseta0->FindBin(fEtaRange[0]));
431 Double_t higheta0 = axiseta0->GetBinUpEdge(axiseta0->FindBin(fEtaRange[1]));
432 Double_t loweta1 = axiseta1->GetBinLowEdge(axiseta1->FindBin(fEtaRange[0]));
433 Double_t higheta1 = axiseta1->GetBinUpEdge(axiseta1->FindBin(fEtaRange[1]));
434 AliInfo(Form("0 Low eta %f and high eta %f\n",loweta0,higheta0));
435 AliInfo(Form("1 Low eta %f and high eta %f\n",loweta1,higheta1));
441 if(charge != kAllCharge) {
442 if((3+((Int_t)(ndimensions/2.))) < ndimensions) {
444 TAxis *axischarge0 = correlationmatrix->GetAxis(3);
445 TAxis *axischarge1 = correlationmatrix->GetAxis(3+((Int_t)(ndimensions/2.)));
447 Int_t bins0 = axischarge0->GetNbins();
448 Int_t bins1 = axischarge1->GetNbins();
450 AliDebug(1, Form("Number of charge bins: %d and %d\n",bins0,bins1));
452 AliError("Problem in the dimensions");
456 axischarge0->SetRangeUser(charge,charge);
457 axischarge1->SetRangeUser(charge,charge);
459 Double_t lowcharge0 = axischarge0->GetBinLowEdge(axischarge0->FindBin(charge));
460 Double_t highcharge0 = axischarge0->GetBinUpEdge(axischarge0->FindBin(charge));
461 Double_t lowcharge1 = axischarge1->GetBinLowEdge(axischarge1->FindBin(charge));
462 Double_t highcharge1 = axischarge1->GetBinUpEdge(axischarge1->FindBin(charge));
463 AliInfo(Form("0 Low charge %f and high charge %f\n",lowcharge0,highcharge0));
464 AliInfo(Form("1 Low charge %f and high charge %f\n",lowcharge1,highcharge1));
470 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
472 Int_t *dim = new Int_t[nDim*2];
473 for(Int_t iter=0; iter < nDim; iter++){
474 dim[iter] = dimensions[iter];
475 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
478 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
484 //___________________________________________________________________________
485 void AliHFECorrectSpectrumBase::CorrectFromTheWidth(TH1D *h1) const {
487 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
490 TAxis *axis = h1->GetXaxis();
491 Int_t nbinX = h1->GetNbinsX();
493 for(Int_t i = 1; i <= nbinX; i++) {
495 Double_t width = axis->GetBinWidth(i);
496 Double_t content = h1->GetBinContent(i);
497 Double_t error = h1->GetBinError(i);
498 h1->SetBinContent(i,content/width);
499 h1->SetBinError(i,error/width);
504 //___________________________________________________________________________
505 void AliHFECorrectSpectrumBase::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
507 // Correct statistical error
510 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
511 Int_t nbinX = h1->GetNbinsX();
513 for(Long_t i = 1; i <= nbinX; i++) {
515 Float_t content = h1->GetBinContent(i);
517 Float_t error = TMath::Sqrt(content);
518 backgroundGrid->SetElementError(bins, error);
522 //_________________________________________________________________________
523 TObject* AliHFECorrectSpectrumBase::GetSpectrum(const AliCFContainer * const c, Int_t step) {
524 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
527 //_________________________________________________________________________
528 TObject* AliHFECorrectSpectrumBase::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
530 // Create efficiency grid and calculate efficiency
536 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
537 eff->CalculateEfficiency(step,step0);