]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEspectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEspectrum.cxx
CommitLineData
c04c80e6 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15//
16// Class for spectrum correction
17// Subtraction of hadronic background, Unfolding of the data and
18// Renormalization done here
19// The following containers have to be set:
20// - Correction framework container for real data
21// - Correction framework container for MC (Efficiency Map)
22// - Correction framework container for background coming from data
23// - Correction framework container for background coming from MC
24//
245877d0 25// Author:
26// Raphaelle Bailhache <R.Bailhache@gsi.de>
27// Markus Fasel <M.Fasel@gsi.de>
c04c80e6 28//
29
30#include <TArrayD.h>
31#include <TH1.h>
32#include <TList.h>
33#include <TObjArray.h>
34#include <TROOT.h>
35#include <TCanvas.h>
36#include <TLegend.h>
37#include <TStyle.h>
38#include <TMath.h>
39#include <TAxis.h>
40#include <TGraphErrors.h>
41#include <TFile.h>
42#include <TPad.h>
67fe7bd0 43#include <TH2D.h>
c2690925 44#include <TF1.h>
c04c80e6 45
3a72645a 46#include "AliPID.h"
c04c80e6 47#include "AliCFContainer.h"
48#include "AliCFDataGrid.h"
49#include "AliCFEffGrid.h"
50#include "AliCFGridSparse.h"
51#include "AliCFUnfolding.h"
52#include "AliLog.h"
53
54#include "AliHFEspectrum.h"
3a72645a 55#include "AliHFEcuts.h"
56#include "AliHFEcontainer.h"
c2690925 57#include "AliHFEtools.h"
c04c80e6 58
59ClassImp(AliHFEspectrum)
60
61//____________________________________________________________
62AliHFEspectrum::AliHFEspectrum(const char *name):
63 TNamed(name, ""),
64 fCFContainers(NULL),
65 fTemporaryObjects(NULL),
66 fCorrelation(NULL),
67 fBackground(NULL),
c2690925 68 fEfficiencyFunction(NULL),
69 fWeightCharm(NULL),
3a72645a 70 fInclusiveSpectrum(kTRUE),
c04c80e6 71 fDumpToFile(kFALSE),
8c1c76e9 72 fEtaSelected(kFALSE),
c2690925 73 fUnSetCorrelatedErrors(kTRUE),
e156c3bb 74 fSetSmoothing(kFALSE),
c2690925 75 fIPanaHadronBgSubtract(kFALSE),
76 fIPanaCharmBgSubtract(kFALSE),
77 fIPanaConversionBgSubtract(kFALSE),
78 fIPanaNonHFEBgSubtract(kFALSE),
79 fNonHFEbgMethod2(kFALSE),
8c1c76e9 80 fNonHFEmode(0),
3a72645a 81 fNbDimensions(1),
c2690925 82 fNMCEvents(0),
83 fNMCbgEvents(0),
c04c80e6 84 fStepMC(-1),
85 fStepTrue(-1),
86 fStepData(-1),
3a72645a 87 fStepBeforeCutsV0(-1),
88 fStepAfterCutsV0(-1),
c04c80e6 89 fStepGuessedUnfolding(-1),
3a72645a 90 fNumberOfIterations(5),
c2690925 91 fChargeChoosen(-1),
92 fNCentralityBinAtTheEnd(0),
93 fHadronEffbyIPcut(NULL),
8c1c76e9 94 fConversionEff(NULL),
95 fNonHFEEff(NULL),
c2690925 96 fBeamType(0),
3a72645a 97 fDebugLevel(0)
c04c80e6 98{
99 //
100 // Default constructor
101 //
c2690925 102
103 for(Int_t k = 0; k < 20; k++){
104 fNEvents[k] = 0;
105 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
106 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
107 }
8c1c76e9 108 memset(fEtaRange, 0, sizeof(Double_t) * 2);
c04c80e6 109}
110
111//____________________________________________________________
112AliHFEspectrum::~AliHFEspectrum(){
113 //
114 // Destructor
115 //
116 if(fCFContainers) delete fCFContainers;
117 if(fTemporaryObjects){
118 fTemporaryObjects->Clear();
119 delete fTemporaryObjects;
120 }
121}
3a72645a 122//____________________________________________________________
c2690925 123Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
3a72645a 124 //
125 // Init what we need for the correction:
126 //
127 // Raw spectrum, hadron contamination
128 // MC efficiency maps, correlation matrix
129 // V0 efficiency if wanted
130 //
131 // This for a given dimension.
132 // If no inclusive spectrum, then take only efficiency map for beauty electron
133 // and the appropriate correlation matrix
134 //
a199006c 135 Int_t kNdim = 3;
c2690925 136 if(fBeamType==0) kNdim=3;
137 if(fBeamType==1) kNdim=4;
138
139 Int_t dims[kNdim];
140 // Get the requested format
141 if(fBeamType==0)
142 {
143 // pp analysis
144 switch(fNbDimensions){
145 case 1: dims[0] = 0;
146 break;
147 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
148 break;
149 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
150 break;
151 default:
152 AliError("Container with this number of dimensions not foreseen (yet)");
153 return kFALSE;
154 };
155 }
156
157 if(fBeamType==1)
158 {
159 // PbPb analysis; centrality as first dimension
160 Int_t nbDimensions = fNbDimensions;
161 fNbDimensions = fNbDimensions + 1;
162 switch(nbDimensions){
163 case 1:
6c70d827 164 dims[0] = 5;
165 dims[1] = 0;
166 break;
c2690925 167 case 2:
6c70d827 168 dims[0] = 5;
169 for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
170 break;
171 case 3:
172 dims[0] = 5;
173 for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
174 break;
c2690925 175 default:
176 AliError("Container with this number of dimensions not foreseen (yet)");
177 return kFALSE;
178 };
179 }
180
181
3a72645a 182
3a72645a 183 // Data container: raw spectrum + hadron contamination
c2690925 184 AliCFContainer *datacontainer = 0x0;
185 if(fInclusiveSpectrum) {
186 datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
187 }
188 else{
189 datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
190 }
3a72645a 191 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
192 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
193
c2690925 194 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
195 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
3a72645a 196 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
c2690925 197
3a72645a 198 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
199 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
200
201 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
202 AliCFContainer *mccontaineresd = 0x0;
203 AliCFContainer *mccontainermc = 0x0;
8c1c76e9 204 AliCFContainer *mccontainermcbg = 0x0;
205 AliCFContainer *nonHFEweightedContainer = 0x0;
206 AliCFContainer *convweightedContainer = 0x0;
207
3a72645a 208 if(fInclusiveSpectrum) {
209 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
210 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
211 }
212 else {
213 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
214 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
8c1c76e9 215 mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
216 nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
217 convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
218 if((!nonHFEweightedContainer) || (!convweightedContainer)) return kFALSE;
3a72645a 219 }
220 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
221
222 Int_t source = -1;
8c1c76e9 223 if(!fInclusiveSpectrum) source = 1; //beauty
c2690925 224 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
225 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
3a72645a 226 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
227 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
228 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
229
c2690925 230 // set charm, nonHFE container to estimate BG
231 if(!fInclusiveSpectrum){
232 source = 0; //charm
8c1c76e9 233 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
c2690925 234 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
8c1c76e9 235
c2690925 236 source = 2; //conversion
8c1c76e9 237 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
238 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
239 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
240 fConversionEff = (TH1D*)efficiencyConv->Project(0);
241
c2690925 242 source = 3; //non HFE except for the conversion
8c1c76e9 243 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
244 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
245 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
246 fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0);
247
248 AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
249 SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
250 AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
251 SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
c2690925 252 }
8c1c76e9 253
254// MC container: correlation matrix
3a72645a 255 THnSparseF *mccorrelation = 0x0;
bf892a6a 256 if(fInclusiveSpectrum) {
257 if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
6c70d827 258 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
259 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
260 else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
261 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
e156c3bb 262
263 if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
bf892a6a 264 }
c2690925 265 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
266 //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
3a72645a 267 if(!mccorrelation) return kFALSE;
268 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
269 if(!mccorrelationD) {
270 printf("No correlation\n");
271 return kFALSE;
272 }
273 SetCorrelation(mccorrelationD);
274
275 // V0 container Electron, pt eta phi
276 if(v0hfecontainer) {
277 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
278 if(!containerV0) return kFALSE;
279 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron+1);
280 if(!containerV0Electron) return kFALSE;
281 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
8c1c76e9 282 }
283
3a72645a 284
285 if(fDebugLevel>0){
286
287 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
288 contaminationspectrum->SetName("contaminationspectrum");
289 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
290 ccontaminationspectrum->Divide(3,1);
291 ccontaminationspectrum->cd(1);
6555e2ad 292 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
293 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
294 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
3a72645a 295 contaminationspectrum2dpteta->SetStats(0);
296 contaminationspectrum2dpteta->SetTitle("");
297 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
298 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
299 contaminationspectrum2dptphi->SetStats(0);
300 contaminationspectrum2dptphi->SetTitle("");
301 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
302 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
303 contaminationspectrum2detaphi->SetStats(0);
304 contaminationspectrum2detaphi->SetTitle("");
305 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
306 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
307 contaminationspectrum2dptphi->Draw("colz");
308 ccontaminationspectrum->cd(2);
309 contaminationspectrum2dpteta->Draw("colz");
310 ccontaminationspectrum->cd(3);
311 contaminationspectrum2detaphi->Draw("colz");
312
c2690925 313 TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
314 ccorrelation->cd(1);
315 if(mccorrelationD) {
316 if(fBeamType==0){
317 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
318 ptcorrelation->Draw("colz");
319 }
320 if(fBeamType==1){
321 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
322 ptcorrelation->Draw("colz");
323 }
324 }
3a72645a 325 }
326
327
328 return kTRUE;
329
330
331}
c04c80e6 332
c2690925 333
c04c80e6 334//____________________________________________________________
3a72645a 335Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
c04c80e6 336 //
337 // Correct the spectrum for efficiency and unfolding
338 // with both method and compare
339 //
340
341 gStyle->SetPalette(1);
342 gStyle->SetOptStat(1111);
343 gStyle->SetPadBorderMode(0);
344 gStyle->SetCanvasColor(10);
345 gStyle->SetPadLeftMargin(0.13);
346 gStyle->SetPadRightMargin(0.13);
67fe7bd0 347
c2690925 348 printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
349
350 ///////////////////////////
351 // Check initialization
352 ///////////////////////////
353
354 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
355 AliInfo("You have to init before");
356 return kFALSE;
357 }
358
a199006c 359 if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
c2690925 360 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
361 return kFALSE;
362 }
363
364 SetNumberOfIteration(10);
365 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
366
367 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
368 //////////////////////////////////
369 // Subtract hadron background
370 /////////////////////////////////
371 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
372 if(subtractcontamination) {
373 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
374 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
375 }
376
377 ////////////////////////////////////////////////
378 // Correct for TPC efficiency from V0
379 ///////////////////////////////////////////////
380 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
381 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
8c1c76e9 382 if(dataContainerV0){printf("Got the V0 container\n");
c2690925 383 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
384 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
385 }
386
387 //////////////////////////////////////////////////////////////////////////////
388 // Correct for efficiency parametrized (if TPC efficiency is parametrized)
389 /////////////////////////////////////////////////////////////////////////////
390 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
391 if(fEfficiencyFunction){
392 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
393 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
394 }
395
396 ///////////////
397 // Unfold
398 //////////////
399 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
400 if(!listunfolded){
401 printf("Unfolded failed\n");
402 return kFALSE;
403 }
404 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
405 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
406 if(!correctedspectrum){
407 AliError("No corrected spectrum\n");
408 return kFALSE;
409 }
410 if(!residualspectrum){
411 AliError("No residul spectrum\n");
412 return kFALSE;
413 }
414
415 /////////////////////
416 // Simply correct
417 ////////////////////
418 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
419
420
421 //////////
422 // Plot
423 //////////
424 if(fDebugLevel > 0.0) {
425
a199006c 426 Int_t ptpr = 0;
c2690925 427 if(fBeamType==0) ptpr=0;
428 if(fBeamType==1) ptpr=1;
429
430 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
431 ccorrected->Divide(2,1);
432 ccorrected->cd(1);
433 gPad->SetLogy();
434 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
435 correctedspectrumD->SetTitle("");
436 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
437 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
438 correctedspectrumD->SetMarkerStyle(26);
439 correctedspectrumD->SetMarkerColor(kBlue);
440 correctedspectrumD->SetLineColor(kBlue);
441 correctedspectrumD->Draw("AP");
442 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
443 alltogetherspectrumD->SetTitle("");
444 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
445 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
446 alltogetherspectrumD->SetMarkerStyle(25);
447 alltogetherspectrumD->SetMarkerColor(kBlack);
448 alltogetherspectrumD->SetLineColor(kBlack);
449 alltogetherspectrumD->Draw("P");
450 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
451 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
452 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
453 legcorrected->Draw("same");
454 ccorrected->cd(2);
455 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
456 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
457 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
458 ratiocorrected->SetName("ratiocorrected");
459 ratiocorrected->SetTitle("");
460 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
461 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
462 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
463 ratiocorrected->SetStats(0);
464 ratiocorrected->Draw();
465
a199006c 466 //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
467 //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
468 //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
469 //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
470
471 TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
472 TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
473 TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
474 TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
c2690925 475
8c1c76e9 476
c2690925 477
478 if(fBeamType==1) {
479
e156c3bb 480 TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
481 ccorrectedallspectra->Divide(2,1);
c2690925 482 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
483 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
484
485 THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
486 TAxis *cenaxisa = sparsesu->GetAxis(0);
487 THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
488 TAxis *cenaxisb = sparsed->GetAxis(0);
489 Int_t nbbin = cenaxisb->GetNbins();
490 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
491 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
492 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
493 TString titlee("corrected_centrality_bin_");
494 titlee += "[";
495 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
496 titlee += "_";
497 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
498 titlee += "[";
499 TString titleec("corrected_check_projection_bin_");
500 titleec += "[";
501 titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
502 titleec += "_";
503 titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
504 titleec += "[";
505 TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
506 titleenameunotnormalized += "[";
507 titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
508 titleenameunotnormalized += "_";
509 titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
510 titleenameunotnormalized += "[";
511 TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
512 titleenameunormalized += "[";
513 titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
514 titleenameunormalized += "_";
515 titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
516 titleenameunormalized += "[";
517 TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
518 titleenamednotnormalized += "[";
519 titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
520 titleenamednotnormalized += "_";
521 titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
522 titleenamednotnormalized += "[";
523 TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
524 titleenamednormalized += "[";
525 titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
526 titleenamednormalized += "_";
527 titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
528 titleenamednormalized += "[";
529 Int_t nbEvents = 0;
530 for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
e156c3bb 531 printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
c2690925 532 nbEvents += fNEvents[k];
533 }
e156c3bb 534 Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
535 Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
536 printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
537 Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
538 Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
539 printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
c2690925 540 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
541 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
542 TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
543 ccorrectedcheck->cd(1);
544 TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
545 TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
546 aftersuc->Draw();
547 aftersdc->Draw("same");
548 TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
549 ccorrectede->Divide(2,1);
550 ccorrectede->cd(1);
551 gPad->SetLogy();
552 TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
553 CorrectFromTheWidth(aftersu);
554 aftersu->SetName((const char*)titleenameunotnormalized);
6c70d827 555 unfoldingspectrac[binc] = *aftersu;
c2690925 556 ccorrectede->cd(1);
557 TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
558 aftersun->SetTitle("");
559 aftersun->GetYaxis()->SetTitleOffset(1.5);
560 aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
561 aftersun->SetMarkerStyle(26);
562 aftersun->SetMarkerColor(kBlue);
563 aftersun->SetLineColor(kBlue);
564 aftersun->Draw("AP");
565 aftersun->SetName((const char*)titleenameunormalized);
6c70d827 566 unfoldingspectracn[binc] = *aftersun;
c2690925 567 ccorrectede->cd(1);
568 TH1D *aftersd = (TH1D *) sparsed->Projection(1);
569 CorrectFromTheWidth(aftersd);
570 aftersd->SetName((const char*)titleenamednotnormalized);
6c70d827 571 correctedspectrac[binc] = *aftersd;
c2690925 572 ccorrectede->cd(1);
573 TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
574 aftersdn->SetTitle("");
575 aftersdn->GetYaxis()->SetTitleOffset(1.5);
576 aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
577 aftersdn->SetMarkerStyle(25);
578 aftersdn->SetMarkerColor(kBlack);
579 aftersdn->SetLineColor(kBlack);
580 aftersdn->Draw("P");
581 aftersdn->SetName((const char*)titleenamednormalized);
6c70d827 582 correctedspectracn[binc] = *aftersdn;
c2690925 583 TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
584 legcorrectedud->AddEntry(aftersun,"Corrected","p");
585 legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
586 legcorrectedud->Draw("same");
e156c3bb 587 ccorrectedallspectra->cd(1);
c2690925 588 gPad->SetLogy();
589 TH1D *aftersunn = (TH1D *) aftersun->Clone();
590 aftersunn->SetMarkerStyle(stylee[binc]);
591 aftersunn->SetMarkerColor(colorr[binc]);
592 if(binc==0) aftersunn->Draw("AP");
593 else aftersunn->Draw("P");
594 legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
e156c3bb 595 ccorrectedallspectra->cd(2);
c2690925 596 gPad->SetLogy();
597 TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
598 aftersdnn->SetMarkerStyle(stylee[binc]);
599 aftersdnn->SetMarkerColor(colorr[binc]);
600 if(binc==0) aftersdnn->Draw("AP");
601 else aftersdnn->Draw("P");
602 legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
603 ccorrectede->cd(2);
604 TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
605 TString titleee("ratiocorrected_bin_");
606 titleee += binc;
607 ratiocorrectedbinc->SetName((const char*) titleee);
608 ratiocorrectedbinc->SetTitle("");
609 ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
610 ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
611 ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
612 ratiocorrectedbinc->SetStats(0);
613 ratiocorrectedbinc->Draw();
614 }
615
e156c3bb 616 ccorrectedallspectra->cd(1);
c2690925 617 legtotal->Draw("same");
e156c3bb 618 ccorrectedallspectra->cd(2);
c2690925 619 legtotalg->Draw("same");
620
621 cenaxisa->SetRange(0,nbbin);
622 cenaxisb->SetRange(0,nbbin);
623
624 }
625
626 // Dump to file if needed
627 if(fDumpToFile) {
628 TFile *out = new TFile("finalSpectrum.root","recreate");
629 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
630 correctedspectrumD->Write();
631 alltogetherspectrumD->SetName("AlltogetherSpectrum");
632 alltogetherspectrumD->Write();
633 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
634 ratiocorrected->Write();
635 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
636 correctedspectrum->Write();
637 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
638 alltogetherCorrection->Write();
639 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
6c70d827 640 unfoldingspectrac[binc].Write();
641 unfoldingspectracn[binc].Write();
642 correctedspectrac[binc].Write();
643 correctedspectracn[binc].Write();
c2690925 644 }
645 out->Close(); delete out;
646 }
647
a199006c 648 if (unfoldingspectrac) delete[] unfoldingspectrac;
649 if (unfoldingspectracn) delete[] unfoldingspectracn;
650 if (correctedspectrac) delete[] correctedspectrac;
651 if (correctedspectracn) delete[] correctedspectracn;
652
017dcb19 653 }
c2690925 654
8c1c76e9 655
656
657
c2690925 658 return kTRUE;
659}
660
661//____________________________________________________________
662Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
663 //
664 // Correct the spectrum for efficiency and unfolding for beauty analysis
665 // with both method and compare
666 //
667
668 gStyle->SetPalette(1);
669 gStyle->SetOptStat(1111);
670 gStyle->SetPadBorderMode(0);
671 gStyle->SetCanvasColor(10);
672 gStyle->SetPadLeftMargin(0.13);
673 gStyle->SetPadRightMargin(0.13);
674
3a72645a 675 ///////////////////////////
676 // Check initialization
677 ///////////////////////////
c04c80e6 678
3a72645a 679 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
680 AliInfo("You have to init before");
681 return kFALSE;
682 }
683
684 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
685 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
686 return kFALSE;
687 }
688
c2690925 689 SetNumberOfIteration(10);
3a72645a 690 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
691
692 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
693 //////////////////////////////////
694 // Subtract hadron background
695 /////////////////////////////////
67fe7bd0 696 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
3a72645a 697 if(subtractcontamination) {
698 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
699 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
700 }
701
c2690925 702 /////////////////////////////////////////////////////////////////////////////////////////
703 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
704 /////////////////////////////////////////////////////////////////////////////////////////
705
8c1c76e9 706 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
3a72645a 707 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
708 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
8c1c76e9 709
710 if(fEfficiencyFunction){
711 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
712 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
713 }
714 else if(dataContainerV0){
3a72645a 715 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
716 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
8c1c76e9 717 }
718
719
3a72645a 720 ///////////////
c04c80e6 721 // Unfold
3a72645a 722 //////////////
723 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
c04c80e6 724 if(!listunfolded){
725 printf("Unfolded failed\n");
3a72645a 726 return kFALSE;
c04c80e6 727 }
728 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
729 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
730 if(!correctedspectrum){
731 AliError("No corrected spectrum\n");
3a72645a 732 return kFALSE;
c04c80e6 733 }
67fe7bd0 734 if(!residualspectrum){
8c1c76e9 735 AliError("No residual spectrum\n");
3a72645a 736 return kFALSE;
c04c80e6 737 }
67fe7bd0 738
3a72645a 739 /////////////////////
c04c80e6 740 // Simply correct
3a72645a 741 ////////////////////
742 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
67fe7bd0 743
3a72645a 744
67fe7bd0 745 //////////
c04c80e6 746 // Plot
747 //////////
3a72645a 748 if(fDebugLevel > 0.0) {
749
750 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
751 ccorrected->Divide(2,1);
752 ccorrected->cd(1);
753 gPad->SetLogy();
754 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
755 correctedspectrumD->SetTitle("");
756 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
757 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
758 correctedspectrumD->SetMarkerStyle(26);
759 correctedspectrumD->SetMarkerColor(kBlue);
760 correctedspectrumD->SetLineColor(kBlue);
761 correctedspectrumD->Draw("AP");
762 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
763 alltogetherspectrumD->SetTitle("");
764 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
765 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
766 alltogetherspectrumD->SetMarkerStyle(25);
767 alltogetherspectrumD->SetMarkerColor(kBlack);
768 alltogetherspectrumD->SetLineColor(kBlack);
769 alltogetherspectrumD->Draw("P");
770 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
771 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
772 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
773 legcorrected->Draw("same");
774 ccorrected->cd(2);
775 TH1D *correctedTH1D = correctedspectrum->Projection(0);
6555e2ad 776 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
3a72645a 777 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
778 ratiocorrected->SetName("ratiocorrected");
779 ratiocorrected->SetTitle("");
780 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
781 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
782 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
783 ratiocorrected->SetStats(0);
784 ratiocorrected->Draw();
785
786
787 // Dump to file if needed
788
789 if(fDumpToFile) {
790
8c1c76e9 791 TFile *out;
792 if(fNonHFEmode == 1){
793 out = new TFile("finalSpectrumLow.root","recreate");
794 }
795 else if(fNonHFEmode == 2){
796 out = new TFile("finalSpectrumUp.root","recreate");
797 }
798 else{
799 out = new TFile("finalSpectrum.root","recreate");
800 }
3a72645a 801 out->cd();
802 //
803 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
804 correctedspectrumD->Write();
805 alltogetherspectrumD->SetName("AlltogetherSpectrum");
806 alltogetherspectrumD->Write();
807 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
808 ratiocorrected->Write();
809 //
810 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
811 correctedspectrum->Write();
812 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
813 alltogetherCorrection->Write();
814 //
8c1c76e9 815 out->Close();
816 delete out;
3a72645a 817 }
818
819
820 }
821
3a72645a 822 return kTRUE;
823}
c2690925 824
3a72645a 825//____________________________________________________________
826AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
827 //
828 // Apply background subtraction
829 //
830
831 // Raw spectrum
832 AliCFContainer *dataContainer = GetContainer(kDataContainer);
833 if(!dataContainer){
834 AliError("Data Container not available");
835 return NULL;
836 }
c2690925 837 printf("Step data: %d\n",fStepData);
3a72645a 838 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
839
840 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
841 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
842
843 // Background Estimate
844 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
845 if(!backgroundContainer){
846 AliError("MC background container not found");
847 return NULL;
848 }
849
c2690925 850 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
3a72645a 851 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
852
c2690925 853 if(!fInclusiveSpectrum){
854 //Background subtraction for IP analysis
8c1c76e9 855 TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0);
856 CorrectFromTheWidth(measuredTH1Draw);
857 TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
858 rawspectra->cd();
859 measuredTH1Draw->SetMarkerStyle(20);
860 measuredTH1Draw->Draw();
c2690925 861 if(fIPanaHadronBgSubtract){
862 // Hadron background
863 printf("Hadron background for IP analysis subtracted!\n");
864 Int_t* bins=new Int_t[1];
865 TH1D* htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
866 bins[0]=htemp->GetNbinsX();
867 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
868 hbgContainer->SetGrid(fHadronEffbyIPcut);
869 backgroundGrid->Multiply(hbgContainer,1);
8c1c76e9 870 // draw raw hadron bg spectra
871 TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0);
872 CorrectFromTheWidth(hadronbg);
873 hadronbg->SetMarkerColor(7);
874 hadronbg->SetMarkerStyle(20);
875 rawspectra->cd();
876 hadronbg->Draw("samep");
877 // subtract hadron contamination
c2690925 878 spectrumSubtracted->Add(backgroundGrid,-1.0);
879 }
880 if(fIPanaCharmBgSubtract){
881 // Charm background
882 printf("Charm background for IP analysis subtracted!\n");
883 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
8c1c76e9 884 // draw charm bg spectra
885 TH1D *charmbg= (TH1D *) charmbgContainer->Project(0);
886 CorrectFromTheWidth(charmbg);
887 charmbg->SetMarkerColor(3);
888 charmbg->SetMarkerStyle(20);
889 rawspectra->cd();
890 charmbg->Draw("samep");
891 // subtract charm background
c2690925 892 spectrumSubtracted->Add(charmbgContainer,-1.0);
893 }
894 if(fIPanaConversionBgSubtract){
895 // Conversion background
896 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
8c1c76e9 897 // draw conversion bg spectra
898 TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0);
899 CorrectFromTheWidth(conversionbg);
900 conversionbg->SetMarkerColor(4);
901 conversionbg->SetMarkerStyle(20);
902 rawspectra->cd();
903 conversionbg->Draw("samep");
904 // subtract conversion background
c2690925 905 spectrumSubtracted->Add(conversionbgContainer,-1.0);
906 printf("Conversion background subtraction is preliminary!\n");
907 }
908 if(fIPanaNonHFEBgSubtract){
909 // NonHFE background
910 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
8c1c76e9 911 // draw Dalitz/dielectron bg spectra
912 TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0);
913 CorrectFromTheWidth(nonhfebg);
914 nonhfebg->SetMarkerColor(6);
915 nonhfebg->SetMarkerStyle(20);
916 rawspectra->cd();
917 nonhfebg->Draw("samep");
918 // subtract Dalitz/dielectron background
c2690925 919 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
920 printf("Non HFE background subtraction is preliminary!\n");
921 }
8c1c76e9 922 TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0);
923 CorrectFromTheWidth(rawbgsubtracted);
924 rawbgsubtracted->SetMarkerStyle(24);
925 rawspectra->cd();
926 rawbgsubtracted->Draw("samep");
c2690925 927 }
928 else{
929 // Subtract
930 spectrumSubtracted->Add(backgroundGrid,-1.0);
931 }
932
3a72645a 933 if(setBackground){
934 if(fBackground) delete fBackground;
935 fBackground = backgroundGrid;
936 } else delete backgroundGrid;
937
938
939 if(fDebugLevel > 0) {
c2690925 940
941 Int_t ptpr;
942 if(fBeamType==0) ptpr=0;
943 if(fBeamType==1) ptpr=1;
67fe7bd0 944
945 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
3a72645a 946 cbackgroundsubtraction->Divide(3,1);
67fe7bd0 947 cbackgroundsubtraction->cd(1);
3a72645a 948 gPad->SetLogy();
c2690925 949 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
950 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
67fe7bd0 951 CorrectFromTheWidth(measuredTH1Daftersubstraction);
952 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
953 measuredTH1Daftersubstraction->SetStats(0);
954 measuredTH1Daftersubstraction->SetTitle("");
955 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
956 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
957 measuredTH1Daftersubstraction->SetMarkerStyle(25);
958 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
959 measuredTH1Daftersubstraction->SetLineColor(kBlack);
960 measuredTH1Dbeforesubstraction->SetStats(0);
961 measuredTH1Dbeforesubstraction->SetTitle("");
962 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
963 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
964 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
965 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
966 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
967 measuredTH1Daftersubstraction->Draw();
968 measuredTH1Dbeforesubstraction->Draw("same");
969 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
3a72645a 970 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
971 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
67fe7bd0 972 legsubstraction->Draw("same");
973 cbackgroundsubtraction->cd(2);
3a72645a 974 gPad->SetLogy();
67fe7bd0 975 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
976 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
977 ratiomeasuredcontamination->SetTitle("");
978 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
979 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
c2690925 980 ratiomeasuredcontamination->Sumw2();
67fe7bd0 981 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
982 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
983 ratiomeasuredcontamination->SetStats(0);
984 ratiomeasuredcontamination->SetMarkerStyle(26);
985 ratiomeasuredcontamination->SetMarkerColor(kBlack);
986 ratiomeasuredcontamination->SetLineColor(kBlack);
c2690925 987 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
988 ratiomeasuredcontamination->SetBinError(k+1,0.0);
989 }
990 ratiomeasuredcontamination->Draw("P");
3a72645a 991 cbackgroundsubtraction->cd(3);
c2690925 992 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
3a72645a 993 CorrectFromTheWidth(measuredTH1background);
994 measuredTH1background->SetStats(0);
995 measuredTH1background->SetTitle("");
996 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
997 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
998 measuredTH1background->SetMarkerStyle(26);
999 measuredTH1background->SetMarkerColor(kBlack);
1000 measuredTH1background->SetLineColor(kBlack);
1001 measuredTH1background->Draw();
c2690925 1002
1003 if(fBeamType==1) {
1004
1005 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1006 cbackgrounde->Divide(2,1);
1007 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1008 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1009
1010 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1011 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1012 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1013 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1014 Int_t nbbin = cenaxisb->GetNbins();
1015 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1016 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1017 for(Int_t binc = 0; binc < nbbin; binc++){
1018 TString titlee("BackgroundSubtraction_centrality_bin_");
1019 titlee += binc;
1020 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1021 cbackground->Divide(2,1);
1022 cbackground->cd(1);
1023 gPad->SetLogy();
1024 cenaxisa->SetRange(binc+1,binc+1);
1025 cenaxisb->SetRange(binc+1,binc+1);
1026 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1027 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1028 CorrectFromTheWidth(aftersubstraction);
1029 CorrectFromTheWidth(beforesubstraction);
1030 aftersubstraction->SetStats(0);
1031 aftersubstraction->SetTitle((const char*)titlee);
1032 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1033 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1034 aftersubstraction->SetMarkerStyle(25);
1035 aftersubstraction->SetMarkerColor(kBlack);
1036 aftersubstraction->SetLineColor(kBlack);
1037 beforesubstraction->SetStats(0);
1038 beforesubstraction->SetTitle((const char*)titlee);
1039 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1040 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1041 beforesubstraction->SetMarkerStyle(24);
1042 beforesubstraction->SetMarkerColor(kBlue);
1043 beforesubstraction->SetLineColor(kBlue);
1044 aftersubstraction->Draw();
1045 beforesubstraction->Draw("same");
1046 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1047 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1048 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1049 lega->Draw("same");
1050 cbackgrounde->cd(1);
1051 gPad->SetLogy();
1052 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1053 aftersubtractionn->SetMarkerStyle(stylee[binc]);
1054 aftersubtractionn->SetMarkerColor(colorr[binc]);
1055 if(binc==0) aftersubtractionn->Draw();
1056 else aftersubtractionn->Draw("same");
1057 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1058 cbackgrounde->cd(2);
1059 gPad->SetLogy();
1060 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1061 aftersubtractionng->SetMarkerStyle(stylee[binc]);
1062 aftersubtractionng->SetMarkerColor(colorr[binc]);
1063 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1064 if(binc==0) aftersubtractionng->Draw();
1065 else aftersubtractionng->Draw("same");
1066 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1067 cbackground->cd(2);
1068 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1069 ratiocontamination->SetName("ratiocontamination");
1070 ratiocontamination->SetTitle((const char*)titlee);
1071 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1072 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1073 ratiocontamination->Add(aftersubstraction,-1.0);
1074 ratiocontamination->Divide(beforesubstraction);
1075 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1076 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1077 ratiocontamination->SetBinError(nbinpt+1,0.0);
1078 }
1079 ratiocontamination->SetStats(0);
1080 ratiocontamination->SetMarkerStyle(26);
1081 ratiocontamination->SetMarkerColor(kBlack);
1082 ratiocontamination->SetLineColor(kBlack);
1083 ratiocontamination->Draw("P");
1084
1085 }
1086
1087 cbackgrounde->cd(1);
1088 legtotal->Draw("same");
1089 cbackgrounde->cd(2);
1090 legtotalg->Draw("same");
1091
1092 cenaxisa->SetRange(0,nbbin);
1093 cenaxisb->SetRange(0,nbbin);
1094
1095 }
1096
3a72645a 1097
67fe7bd0 1098 }
1099
3a72645a 1100 return spectrumSubtracted;
c04c80e6 1101}
c2690925 1102
1103//____________________________________________________________
1104AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1105 //
1106 // calculate charm background
1107 //
1108
1109 Double_t evtnorm=0;
8c1c76e9 1110 if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents);
c2690925 1111
1112 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1113 if(!mcContainer){
1114 AliError("MC Container not available");
1115 return NULL;
1116 }
1117
1118 if(!fCorrelation){
1119 AliError("No Correlation map available");
1120 return NULL;
1121 }
1122
c2690925 1123 AliCFDataGrid *charmBackgroundGrid= 0x0;
8c1c76e9 1124 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1125 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1126
c2690925 1127 Int_t* bins=new Int_t[1];
8c1c76e9 1128 bins[0]=charmbgaftertofpid->GetNbinsX();
1129 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins);
1130 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1131 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
1132
1133 charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
1134 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
1135
c2690925 1136 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
8c1c76e9 1137 weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1138 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
1139 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1140 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
c2690925 1141
1142 // Efficiency (set efficiency to 1 for only folding)
1143 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1144 efficiencyD->CalculateEfficiency(0,0);
1145
1146 // Folding
8c1c76e9 1147 Int_t nDim = 1;
c2690925 1148 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1149 folding.SetMaxNumberOfIterations(1);
1150 folding.Unfold();
1151
1152 // Results
1153 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1154 THnSparse* result=(THnSparse*)result1->Clone();
1155 charmBackgroundGrid->SetGrid(result);
8c1c76e9 1156 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
1157
1158 //Charm background evaluation plots
1159
1160 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1161 cCharmBgEval->Divide(3,1);
1162
1163 cCharmBgEval->cd(1);
1164 charmbgaftertofpid->Scale(evtnorm);
1165 CorrectFromTheWidth(charmbgaftertofpid);
1166 charmbgaftertofpid->SetMarkerStyle(25);
1167 charmbgaftertofpid->Draw("p");
1168 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1169 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1170 gPad->SetLogy();
1171
1172 CorrectFromTheWidth(charmbgafteripcut);
1173 charmbgafteripcut->SetMarkerStyle(24);
1174 charmbgafteripcut->Draw("samep");
1175
1176 CorrectFromTheWidth(charmbgafterweight);
1177 charmbgafterweight->SetMarkerStyle(24);
1178 charmbgafterweight->SetMarkerColor(4);
1179 charmbgafterweight->Draw("samep");
1180
1181 CorrectFromTheWidth(charmbgafterfolding);
1182 charmbgafterfolding->SetMarkerStyle(24);
1183 charmbgafterfolding->SetMarkerColor(2);
1184 charmbgafterfolding->Draw("samep");
1185
1186 cCharmBgEval->cd(2);
1187 parametrizedcharmpidipeff->SetMarkerStyle(24);
1188 parametrizedcharmpidipeff->Draw("p");
1189 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1190
1191 cCharmBgEval->cd(3);
1192 charmweightingfc->SetMarkerStyle(24);
1193 charmweightingfc->Draw("p");
1194 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1195 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1196
1197 cCharmBgEval->cd(1);
1198 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1199 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1200 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1201 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1202 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1203 legcharmbg->Draw("same");
1204
1205 cCharmBgEval->cd(2);
1206 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1207 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1208 legcharmbg2->Draw("same");
1209
1210 CorrectStatErr(charmBackgroundGrid);
c2690925 1211
1212 return charmBackgroundGrid;
1213}
1214
1215//____________________________________________________________
1216AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1217 //
1218 // calculate conversion background
1219 //
1220
a199006c 1221 Double_t evtnorm[1] = {0.0};
c2690925 1222 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1223 printf("check event!!! %lf \n",evtnorm[0]);
1224
1225 // Background Estimate
8c1c76e9 1226 AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
c2690925 1227 if(!backgroundContainer){
1228 AliError("MC background container not found");
1229 return NULL;
1230 }
1231
8c1c76e9 1232 Int_t stepbackground = 3;
c2690925 1233 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1234 backgroundGrid->Scale(evtnorm);
1235
8c1c76e9 1236 Int_t* bins=new Int_t[1];
1237 bins[0]=fConversionEff->GetNbinsX();
1238 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1239 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1240 backgroundGrid->Multiply(weightedConversionContainer,1.0);
1241
1242 CorrectStatErr(backgroundGrid);
1243
c2690925 1244 return backgroundGrid;
1245}
1246
1247
1248//____________________________________________________________
1249AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1250 //
8c1c76e9 1251 // calculate non-HFE background
c2690925 1252 //
1253
a199006c 1254 Double_t evtnorm[1] = {0.0};
8c1c76e9 1255 if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1256 printf("check event!!! %lf \n",evtnorm[0]);
c2690925 1257
8c1c76e9 1258 // Background Estimate
1259 AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1260 if(!backgroundContainer){
1261 AliError("MC background container not found");
1262 return NULL;
c2690925 1263 }
c2690925 1264
c2690925 1265
8c1c76e9 1266 Int_t stepbackground = 3;
1267 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1268 backgroundGrid->Scale(evtnorm);
c2690925 1269
8c1c76e9 1270 Int_t* bins=new Int_t[1];
1271 bins[0]=fNonHFEEff->GetNbinsX();
1272 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1273 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1274 backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
c2690925 1275
8c1c76e9 1276 CorrectStatErr(backgroundGrid);
c2690925 1277
8c1c76e9 1278 return backgroundGrid;
c2690925 1279}
1280
1281//____________________________________________________________
1282AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1283
1284 //
1285 // Apply TPC pid efficiency correction from parametrisation
1286 //
1287
1288 // Data in the right format
1289 AliCFDataGrid *dataGrid = 0x0;
1290 if(bgsubpectrum) {
1291 dataGrid = bgsubpectrum;
1292 }
1293 else {
1294
1295 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1296 if(!dataContainer){
1297 AliError("Data Container not available");
1298 return NULL;
1299 }
1300
1301 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1302 }
1303
1304 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1305 result->SetName("ParametrizedEfficiencyBefore");
1306 THnSparse *h = result->GetGrid();
1307 Int_t nbdimensions = h->GetNdimensions();
1308 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1309
1310 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1311 if(!dataContainer){
1312 AliError("Data Container not available");
1313 return NULL;
1314 }
1315 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1316 dataContainerbis->Add(dataContainerbis,-1.0);
1317
1318
1319 Int_t* coord = new Int_t[nbdimensions];
1320 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1321 Double_t* points = new Double_t[nbdimensions];
1322
1323
1324 ULong64_t nEntries = h->GetNbins();
1325 for (ULong64_t i = 0; i < nEntries; ++i) {
1326
1327 Double_t value = h->GetBinContent(i, coord);
1328 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1329 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1330
1331 // Get the bin co-ordinates given an coord
1332 for (Int_t j = 0; j < nbdimensions; ++j)
1333 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1334
1335 if (!fEfficiencyFunction->IsInside(points))
1336 continue;
1337 TF1::RejectPoint(kFALSE);
1338
1339 // Evaulate function at points
1340 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1341 //printf("Value efficiency is %f\n",valueEfficiency);
1342
1343 if(valueEfficiency > 0.0) {
1344 h->SetBinContent(coord,value/valueEfficiency);
1345 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1346 }
1347 Double_t error = h->GetBinError(i);
1348 h->SetBinError(coord,error/valueEfficiency);
1349 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1350
1351
1352 }
1353
6c70d827 1354 delete[] coord;
1355 delete[] points;
1356
c2690925 1357 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1358
1359 if(fDebugLevel > 0) {
1360
1361 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1362 cEfficiencyParametrized->Divide(2,1);
1363 cEfficiencyParametrized->cd(1);
1364 TH1D *afterE = (TH1D *) resultt->Project(0);
1365 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1366 CorrectFromTheWidth(afterE);
1367 CorrectFromTheWidth(beforeE);
1368 afterE->SetStats(0);
1369 afterE->SetTitle("");
1370 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1371 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1372 afterE->SetMarkerStyle(25);
1373 afterE->SetMarkerColor(kBlack);
1374 afterE->SetLineColor(kBlack);
1375 beforeE->SetStats(0);
1376 beforeE->SetTitle("");
1377 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1378 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1379 beforeE->SetMarkerStyle(24);
1380 beforeE->SetMarkerColor(kBlue);
1381 beforeE->SetLineColor(kBlue);
1382 gPad->SetLogy();
1383 afterE->Draw();
1384 beforeE->Draw("same");
1385 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1386 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1387 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1388 legefficiencyparametrized->Draw("same");
1389 cEfficiencyParametrized->cd(2);
1390 fEfficiencyFunction->Draw();
1391 //cEfficiencyParametrized->cd(3);
1392 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1393 //ratioefficiency->Divide(afterE);
1394 //ratioefficiency->Draw();
1395
1396
1397 }
1398
1399
1400 return resultt;
1401
1402}
c04c80e6 1403//____________________________________________________________
3a72645a 1404AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1405
c04c80e6 1406 //
3a72645a 1407 // Apply TPC pid efficiency correction from V0
c04c80e6 1408 //
1409
3a72645a 1410 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1411 if(!v0Container){
1412 AliError("V0 Container not available");
c04c80e6 1413 return NULL;
1414 }
3a72645a 1415
1416 // Efficiency
1417 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1418 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1419
1420 // Data in the right format
1421 AliCFDataGrid *dataGrid = 0x0;
1422 if(bgsubpectrum) {
1423 dataGrid = bgsubpectrum;
c04c80e6 1424 }
3a72645a 1425 else {
1426
1427 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1428 if(!dataContainer){
1429 AliError("Data Container not available");
1430 return NULL;
1431 }
c04c80e6 1432
3a72645a 1433 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1434 }
c04c80e6 1435
3a72645a 1436 // Correct
1437 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1438 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1439
3a72645a 1440 if(fDebugLevel > 0) {
c2690925 1441
1442 Int_t ptpr;
1443 if(fBeamType==0) ptpr=0;
1444 if(fBeamType==1) ptpr=1;
3a72645a 1445
1446 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1447 cV0Efficiency->Divide(2,1);
1448 cV0Efficiency->cd(1);
c2690925 1449 TH1D *afterE = (TH1D *) result->Project(ptpr);
1450 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
3a72645a 1451 afterE->SetStats(0);
1452 afterE->SetTitle("");
1453 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1454 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1455 afterE->SetMarkerStyle(25);
1456 afterE->SetMarkerColor(kBlack);
1457 afterE->SetLineColor(kBlack);
1458 beforeE->SetStats(0);
1459 beforeE->SetTitle("");
1460 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1461 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1462 beforeE->SetMarkerStyle(24);
1463 beforeE->SetMarkerColor(kBlue);
1464 beforeE->SetLineColor(kBlue);
1465 afterE->Draw();
1466 beforeE->Draw("same");
1467 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1468 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1469 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1470 legV0efficiency->Draw("same");
1471 cV0Efficiency->cd(2);
c2690925 1472 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
3a72645a 1473 efficiencyDproj->SetTitle("");
1474 efficiencyDproj->SetStats(0);
1475 efficiencyDproj->SetMarkerStyle(25);
1476 efficiencyDproj->Draw();
c04c80e6 1477
c2690925 1478 if(fBeamType==1) {
1479
1480 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1481 cV0Efficiencye->Divide(2,1);
1482 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1483 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1484
1485 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1486 TAxis *cenaxisa = sparseafter->GetAxis(0);
1487 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1488 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1489 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1490 TAxis *cenaxisc = efficiencya->GetAxis(0);
1491 Int_t nbbin = cenaxisb->GetNbins();
1492 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1493 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1494 for(Int_t binc = 0; binc < nbbin; binc++){
1495 TString titlee("V0Efficiency_centrality_bin_");
1496 titlee += binc;
1497 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1498 ccV0Efficiency->Divide(2,1);
1499 ccV0Efficiency->cd(1);
1500 gPad->SetLogy();
1501 cenaxisa->SetRange(binc+1,binc+1);
1502 cenaxisb->SetRange(binc+1,binc+1);
1503 cenaxisc->SetRange(binc+1,binc+1);
1504 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1505 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1506 CorrectFromTheWidth(aftere);
1507 CorrectFromTheWidth(beforee);
1508 aftere->SetStats(0);
1509 aftere->SetTitle((const char*)titlee);
1510 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1511 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1512 aftere->SetMarkerStyle(25);
1513 aftere->SetMarkerColor(kBlack);
1514 aftere->SetLineColor(kBlack);
1515 beforee->SetStats(0);
1516 beforee->SetTitle((const char*)titlee);
1517 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1518 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1519 beforee->SetMarkerStyle(24);
1520 beforee->SetMarkerColor(kBlue);
1521 beforee->SetLineColor(kBlue);
1522 aftere->Draw();
1523 beforee->Draw("same");
1524 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1525 lega->AddEntry(beforee,"Before correction","p");
1526 lega->AddEntry(aftere,"After correction","p");
1527 lega->Draw("same");
1528 cV0Efficiencye->cd(1);
1529 gPad->SetLogy();
1530 TH1D *afteree = (TH1D *) aftere->Clone();
1531 afteree->SetMarkerStyle(stylee[binc]);
1532 afteree->SetMarkerColor(colorr[binc]);
1533 if(binc==0) afteree->Draw();
1534 else afteree->Draw("same");
1535 legtotal->AddEntry(afteree,(const char*) titlee,"p");
1536 cV0Efficiencye->cd(2);
1537 gPad->SetLogy();
1538 TH1D *aftereeu = (TH1D *) aftere->Clone();
1539 aftereeu->SetMarkerStyle(stylee[binc]);
1540 aftereeu->SetMarkerColor(colorr[binc]);
1541 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1542 if(binc==0) aftereeu->Draw();
1543 else aftereeu->Draw("same");
1544 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1545 ccV0Efficiency->cd(2);
1546 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1547 efficiencyDDproj->SetTitle("");
1548 efficiencyDDproj->SetStats(0);
1549 efficiencyDDproj->SetMarkerStyle(25);
1550 efficiencyDDproj->Draw();
1551
1552 }
1553
1554 cV0Efficiencye->cd(1);
1555 legtotal->Draw("same");
1556 cV0Efficiencye->cd(2);
1557 legtotalg->Draw("same");
1558
1559 cenaxisa->SetRange(0,nbbin);
1560 cenaxisb->SetRange(0,nbbin);
1561 cenaxisc->SetRange(0,nbbin);
1562
1563 }
3a72645a 1564
1565 }
1566
1567
1568 return result;
1569
1570}
c04c80e6 1571//____________________________________________________________
3a72645a 1572TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1573
1574 //
1575 // Unfold and eventually correct for efficiency the bgsubspectrum
1576 //
1577
3a72645a 1578 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
c04c80e6 1579 if(!mcContainer){
1580 AliError("MC Container not available");
1581 return NULL;
1582 }
1583
1584 if(!fCorrelation){
1585 AliError("No Correlation map available");
1586 return NULL;
1587 }
1588
3a72645a 1589 // Data
c04c80e6 1590 AliCFDataGrid *dataGrid = 0x0;
1591 if(bgsubpectrum) {
c04c80e6 1592 dataGrid = bgsubpectrum;
c04c80e6 1593 }
1594 else {
1595
1596 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1597 if(!dataContainer){
1598 AliError("Data Container not available");
1599 return NULL;
1600 }
1601
3a72645a 1602 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1603 }
1604
c04c80e6 1605 // Guessed
3a72645a 1606 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
c04c80e6 1607 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1608
1609 // Efficiency
3a72645a 1610 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1611 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
8c1c76e9 1612
1613 // Consider parameterized IP cut efficiency
1614 if(!fInclusiveSpectrum){
1615 Int_t* bins=new Int_t[1];
1616 bins[0]=44;
1617
1618 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1619 beffContainer->SetGrid(GetBeautyIPEff());
1620 efficiencyD->Multiply(beffContainer,1);
1621 }
1622
c04c80e6 1623
1624 // Unfold
1625
3a72645a 1626 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
c2690925 1627 //unfolding.SetUseCorrelatedErrors();
1628 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
c04c80e6 1629 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
e156c3bb 1630 if(fSetSmoothing) unfolding.UseSmoothing();
c04c80e6 1631 unfolding.Unfold();
1632
1633 // Results
1634 THnSparse* result = unfolding.GetUnfolded();
1635 THnSparse* residual = unfolding.GetEstMeasured();
1636 TList *listofresults = new TList;
1637 listofresults->SetOwner();
1638 listofresults->AddAt((THnSparse*)result->Clone(),0);
1639 listofresults->AddAt((THnSparse*)residual->Clone(),1);
1640
3a72645a 1641 if(fDebugLevel > 0) {
c2690925 1642
1643 Int_t ptpr;
1644 if(fBeamType==0) ptpr=0;
1645 if(fBeamType==1) ptpr=1;
3a72645a 1646
1647 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1648 cresidual->Divide(2,1);
1649 cresidual->cd(1);
1650 gPad->SetLogy();
1651 TGraphErrors* residualspectrumD = Normalize(residual);
1652 if(!residualspectrumD) {
1653 AliError("Number of Events not set for the normalization");
1654 return kFALSE;
1655 }
1656 residualspectrumD->SetTitle("");
1657 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
1658 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1659 residualspectrumD->SetMarkerStyle(26);
1660 residualspectrumD->SetMarkerColor(kBlue);
1661 residualspectrumD->SetLineColor(kBlue);
1662 residualspectrumD->Draw("AP");
1663 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
1664 dataGridBis->SetName("dataGridBis");
1665 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
1666 measuredspectrumD->SetTitle("");
1667 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
1668 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1669 measuredspectrumD->SetMarkerStyle(25);
1670 measuredspectrumD->SetMarkerColor(kBlack);
1671 measuredspectrumD->SetLineColor(kBlack);
1672 measuredspectrumD->Draw("P");
1673 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
1674 legres->AddEntry(residualspectrumD,"Residual","p");
1675 legres->AddEntry(measuredspectrumD,"Measured","p");
1676 legres->Draw("same");
1677 cresidual->cd(2);
c2690925 1678 TH1D *residualTH1D = residual->Projection(ptpr);
1679 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
3a72645a 1680 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
1681 ratioresidual->SetName("ratioresidual");
1682 ratioresidual->SetTitle("");
1683 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
1684 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1685 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
1686 ratioresidual->SetStats(0);
1687 ratioresidual->Draw();
1688
1689 }
1690
c04c80e6 1691 return listofresults;
1692
1693}
1694
1695//____________________________________________________________
3a72645a 1696AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
c04c80e6 1697
1698 //
1699 // Apply unfolding and efficiency correction together to bgsubspectrum
1700 //
1701
3a72645a 1702 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
c04c80e6 1703 if(!mcContainer){
1704 AliError("MC Container not available");
1705 return NULL;
1706 }
1707
c04c80e6 1708 // Efficiency
3a72645a 1709 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
c04c80e6 1710 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1711
8c1c76e9 1712 // Consider parameterized IP cut efficiency
1713 if(!fInclusiveSpectrum){
1714 Int_t* bins=new Int_t[1];
1715 bins[0]=44;
1716
1717 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1718 beffContainer->SetGrid(GetBeautyIPEff());
1719 efficiencyD->Multiply(beffContainer,1);
1720 }
1721
c04c80e6 1722 // Data in the right format
1723 AliCFDataGrid *dataGrid = 0x0;
1724 if(bgsubpectrum) {
c04c80e6 1725 dataGrid = bgsubpectrum;
c04c80e6 1726 }
1727 else {
3a72645a 1728
c04c80e6 1729 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1730 if(!dataContainer){
1731 AliError("Data Container not available");
1732 return NULL;
1733 }
1734
3a72645a 1735 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
c04c80e6 1736 }
1737
1738 // Correct
1739 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1740 result->ApplyEffCorrection(*efficiencyD);
c04c80e6 1741
3a72645a 1742 if(fDebugLevel > 0) {
c2690925 1743
1744 Int_t ptpr;
1745 if(fBeamType==0) ptpr=0;
1746 if(fBeamType==1) ptpr=1;
3a72645a 1747
bf892a6a 1748 printf("Step MC: %d\n",fStepMC);
1749 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1750 printf("Step MC true: %d\n",fStepTrue);
3a72645a 1751 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1752 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
1753 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
1754
1755 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
1756
1757 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
1758 cefficiency->cd(1);
c2690925 1759 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
3a72645a 1760 efficiencymcPIDD->SetTitle("");
1761 efficiencymcPIDD->SetStats(0);
1762 efficiencymcPIDD->SetMarkerStyle(25);
1763 efficiencymcPIDD->Draw();
c2690925 1764 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
3a72645a 1765 efficiencymctrackinggeoD->SetTitle("");
1766 efficiencymctrackinggeoD->SetStats(0);
1767 efficiencymctrackinggeoD->SetMarkerStyle(26);
1768 efficiencymctrackinggeoD->Draw("same");
c2690925 1769 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
3a72645a 1770 efficiencymcallD->SetTitle("");
1771 efficiencymcallD->SetStats(0);
1772 efficiencymcallD->SetMarkerStyle(27);
1773 efficiencymcallD->Draw("same");
c2690925 1774 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
3a72645a 1775 efficiencyesdallD->SetTitle("");
1776 efficiencyesdallD->SetStats(0);
1777 efficiencyesdallD->SetMarkerStyle(24);
1778 efficiencyesdallD->Draw("same");
1779 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
1780 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
1781 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
1782 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
1783 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
1784 legeff->Draw("same");
c2690925 1785
1786 if(fBeamType==1) {
1787
1788 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1789 TAxis *cenaxisa = sparseafter->GetAxis(0);
1790 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1791 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1792 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1793 TAxis *cenaxisc = efficiencya->GetAxis(0);
1794 //Int_t nbbin = cenaxisb->GetNbins();
1795 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1796 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1797 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
1798 TString titlee("Efficiency_centrality_bin_");
1799 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
1800 titlee += "_";
1801 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
1802 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1803 cefficiencye->Divide(2,1);
1804 cefficiencye->cd(1);
1805 gPad->SetLogy();
1806 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1807 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1808 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
1809 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
1810 CorrectFromTheWidth(afterefficiency);
1811 CorrectFromTheWidth(beforeefficiency);
1812 afterefficiency->SetStats(0);
1813 afterefficiency->SetTitle((const char*)titlee);
1814 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1815 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1816 afterefficiency->SetMarkerStyle(25);
1817 afterefficiency->SetMarkerColor(kBlack);
1818 afterefficiency->SetLineColor(kBlack);
1819 beforeefficiency->SetStats(0);
1820 beforeefficiency->SetTitle((const char*)titlee);
1821 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1822 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1823 beforeefficiency->SetMarkerStyle(24);
1824 beforeefficiency->SetMarkerColor(kBlue);
1825 beforeefficiency->SetLineColor(kBlue);
1826 afterefficiency->Draw();
1827 beforeefficiency->Draw("same");
1828 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1829 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
1830 lega->AddEntry(afterefficiency,"After efficiency correction","p");
1831 lega->Draw("same");
1832 cefficiencye->cd(2);
1833 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
1834 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
1835 efficiencyDDproj->SetTitle("");
1836 efficiencyDDproj->SetStats(0);
1837 efficiencyDDproj->SetMarkerStyle(25);
1838 efficiencyDDproj->SetMarkerColor(2);
1839 efficiencyDDproj->Draw();
1840 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
1841 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
1842 efficiencyDDproja->SetTitle("");
1843 efficiencyDDproja->SetStats(0);
1844 efficiencyDDproja->SetMarkerStyle(26);
1845 efficiencyDDproja->SetMarkerColor(4);
1846 efficiencyDDproja->Draw("same");
1847
1848 }
1849 }
1850
1851
3a72645a 1852 }
1853
c04c80e6 1854 return result;
1855
1856}
3a72645a 1857
c04c80e6 1858//__________________________________________________________________________________
c2690925 1859TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
c04c80e6 1860 //
1861 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1862 // Give the final pt spectrum to be compared
1863 //
1864
c2690925 1865 if(fNEvents[i] > 0) {
1866
a199006c 1867 Int_t ptpr = 0;
c2690925 1868 if(fBeamType==0) ptpr=0;
1869 if(fBeamType==1) ptpr=1;
c04c80e6 1870
c2690925 1871 TH1D* projection = spectrum->Projection(ptpr);
c04c80e6 1872 CorrectFromTheWidth(projection);
c2690925 1873 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1874 return graphError;
1875
1876 }
1877
1878 return 0x0;
1879
1880
1881}
1882//__________________________________________________________________________________
c2690925 1883TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
c04c80e6 1884 //
1885 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1886 // Give the final pt spectrum to be compared
1887 //
c2690925 1888 if(fNEvents[i] > 0) {
c04c80e6 1889
a199006c 1890 Int_t ptpr=0;
c2690925 1891 if(fBeamType==0) ptpr=0;
1892 if(fBeamType==1) ptpr=1;
1893
1894 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
c04c80e6 1895 CorrectFromTheWidth(projection);
c2690925 1896 TGraphErrors *graphError = NormalizeTH1(projection,i);
c04c80e6 1897 return graphError;
1898
1899 }
1900
1901 return 0x0;
1902
1903
1904}
1905//__________________________________________________________________________________
c2690925 1906TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
1907 //
1908 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1909 // Give the final pt spectrum to be compared
1910 //
8c1c76e9 1911 Double_t chargecoefficient = 0.5;
1912 if(fChargeChoosen >= 0) chargecoefficient = 1.0;
1913
1914 Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
1915 printf("Normalizing Eta Range %f\n", etarange);
c2690925 1916 if(fNEvents[i] > 0) {
1917
1918 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1919 Double_t p = 0, dp = 0; Int_t point = 1;
1920 Double_t n = 0, dN = 0;
1921 Double_t nCorr = 0, dNcorr = 0;
1922 Double_t errdN = 0, errdp = 0;
1923 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1924 point = ibin - input->GetXaxis()->GetFirst();
1925 p = input->GetXaxis()->GetBinCenter(ibin);
1926 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1927 n = input->GetBinContent(ibin);
1928 dN = input->GetBinError(ibin);
1929
1930 // New point
8c1c76e9 1931 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
c2690925 1932 errdN = 1./(2. * TMath::Pi() * p);
1933 errdp = 1./(2. * TMath::Pi() * p*p) * n;
8c1c76e9 1934 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c2690925 1935
1936 spectrumNormalized->SetPoint(point, p, nCorr);
1937 spectrumNormalized->SetPointError(point, dp, dNcorr);
1938 }
1939 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1940 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1941 spectrumNormalized->SetMarkerStyle(22);
1942 spectrumNormalized->SetMarkerColor(kBlue);
1943 spectrumNormalized->SetLineColor(kBlue);
1944
1945 return spectrumNormalized;
1946
1947 }
1948
1949 return 0x0;
1950
1951
1952}
1953//__________________________________________________________________________________
1954TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
c04c80e6 1955 //
1956 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1957 // Give the final pt spectrum to be compared
1958 //
8c1c76e9 1959 Double_t chargecoefficient = 0.5;
1960 if(fChargeChoosen >= 0) chargecoefficient = 1.0;
1961
1962 Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
1963 printf("Normalizing Eta Range %f\n", etarange);
c2690925 1964 if(normalization > 0) {
c04c80e6 1965
1966 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1967 Double_t p = 0, dp = 0; Int_t point = 1;
1968 Double_t n = 0, dN = 0;
1969 Double_t nCorr = 0, dNcorr = 0;
1970 Double_t errdN = 0, errdp = 0;
1971 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1972 point = ibin - input->GetXaxis()->GetFirst();
1973 p = input->GetXaxis()->GetBinCenter(ibin);
1974 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1975 n = input->GetBinContent(ibin);
1976 dN = input->GetBinError(ibin);
1977
1978 // New point
8c1c76e9 1979 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
c04c80e6 1980 errdN = 1./(2. * TMath::Pi() * p);
1981 errdp = 1./(2. * TMath::Pi() * p*p) * n;
8c1c76e9 1982 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
c04c80e6 1983
1984 spectrumNormalized->SetPoint(point, p, nCorr);
1985 spectrumNormalized->SetPointError(point, dp, dNcorr);
1986 }
1987 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1988 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1989 spectrumNormalized->SetMarkerStyle(22);
1990 spectrumNormalized->SetMarkerColor(kBlue);
1991 spectrumNormalized->SetLineColor(kBlue);
1992
1993 return spectrumNormalized;
1994
1995 }
1996
1997 return 0x0;
1998
1999
2000}
2001//____________________________________________________________
2002void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2003 //
2004 // Set the container for a given step to the
2005 //
2006 if(!fCFContainers) fCFContainers = new TList;
2007 fCFContainers->AddAt(cont, type);
2008}
2009
2010//____________________________________________________________
2011AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2012 //
2013 // Get Correction Framework Container for given type
2014 //
2015 if(!fCFContainers) return NULL;
2016 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2017}
c04c80e6 2018//____________________________________________________________
c2690925 2019AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Int_t positivenegative) {
c04c80e6 2020 //
3a72645a 2021 // Slice bin for a given source of electron
c2690925 2022 // nDim is the number of dimension the corrections are done
2023 // dimensions are the definition of the dimensions
2024 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2025 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
c04c80e6 2026 //
2027
2028 Double_t *varMin = new Double_t[container->GetNVar()],
2029 *varMax = new Double_t[container->GetNVar()];
2030
2031 Double_t *binLimits;
2032 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2033
2034 binLimits = new Double_t[container->GetNBins(ivar)+1];
2035 container->GetBinLimits(ivar,binLimits);
c2690925 2036 varMin[ivar] = binLimits[0];
2037 varMax[ivar] = binLimits[container->GetNBins(ivar)];
2038 // source
2039 if(ivar == 4){
2040 if((source>= 0) && (source<container->GetNBins(ivar))) {
2041 varMin[ivar] = binLimits[source];
2042 varMax[ivar] = binLimits[source];
2043 }
3a72645a 2044 }
c2690925 2045 // charge
2046 if(ivar == 3) {
2047 if((positivenegative>= 0) && (positivenegative<container->GetNBins(ivar))) {
2048 varMin[ivar] = binLimits[positivenegative];
2049 varMax[ivar] = binLimits[positivenegative];
2050 }
3a72645a 2051 }
8c1c76e9 2052 // eta
2053 if(ivar == 1){
2054 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
2055 if(fEtaSelected){
2056 varMin[ivar] = fEtaRange[0];
2057 varMax[ivar] = fEtaRange[1];
2058 }
2059 }
c2690925 2060
3a72645a 2061
c04c80e6 2062 delete[] binLimits;
2063 }
2064
2065 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2066 AddTemporaryObject(k);
2067 delete[] varMin; delete[] varMax;
2068
2069 return k;
2070
2071}
2072
2073//_________________________________________________________________________
3a72645a 2074THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
c04c80e6 2075 //
3a72645a 2076 // Slice correlation
c04c80e6 2077 //
2078
3a72645a 2079 Int_t ndimensions = correlationmatrix->GetNdimensions();
c2690925 2080 //printf("Number of dimension %d correlation map\n",ndimensions);
c04c80e6 2081 if(ndimensions < (2*nDim)) {
2082 AliError("Problem in the dimensions");
2083 return NULL;
2084 }
2085 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
c2690925 2086 //printf("Number of dimension %d container\n",ndimensionsContainer);
c04c80e6 2087
2088 Int_t *dim = new Int_t[nDim*2];
2089 for(Int_t iter=0; iter < nDim; iter++){
2090 dim[iter] = dimensions[iter];
2091 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
c2690925 2092 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
c04c80e6 2093 }
2094
3a72645a 2095 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
c04c80e6 2096
2097 delete[] dim;
2098 return k;
2099
2100}
2101//___________________________________________________________________________
2102void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2103 //
2104 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2105 //
2106
2107 TAxis *axis = h1->GetXaxis();
2108 Int_t nbinX = h1->GetNbinsX();
2109
2110 for(Int_t i = 1; i <= nbinX; i++) {
2111
2112 Double_t width = axis->GetBinWidth(i);
2113 Double_t content = h1->GetBinContent(i);
2114 Double_t error = h1->GetBinError(i);
2115 h1->SetBinContent(i,content/width);
2116 h1->SetBinError(i,error/width);
2117 }
2118
2119}
8c1c76e9 2120
2121//___________________________________________________________________________
2122void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2123 //
2124 // Correct statistical error
2125 //
2126
2127 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2128 Int_t nbinX = h1->GetNbinsX();
2129 Int_t bins[1];
2130 for(Long_t i = 1; i <= nbinX; i++) {
2131 bins[0] = i;
2132 Float_t content = h1->GetBinContent(i);
2133 if(content>0){
2134 Float_t error = TMath::Sqrt(content);
2135 backgroundGrid->SetElementError(bins, error);
2136 }
2137 }
2138}
2139
c04c80e6 2140//____________________________________________________________
2141void AliHFEspectrum::AddTemporaryObject(TObject *o){
2142 //
2143 // Emulate garbage collection on class level
2144 //
2145 if(!fTemporaryObjects) {
2146 fTemporaryObjects= new TList;
2147 fTemporaryObjects->SetOwner();
2148 }
2149 fTemporaryObjects->Add(o);
2150}
2151
2152//____________________________________________________________
2153void AliHFEspectrum::ClearObject(TObject *o){
2154 //
2155 // Do a safe deletion
2156 //
2157 if(fTemporaryObjects){
2158 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2159 delete o;
2160 } else{
2161 // just delete
2162 delete o;
2163 }
2164}
2165//_________________________________________________________________________
c2690925 2166TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
245877d0 2167 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
c04c80e6 2168 return data;
2169}
2170//_________________________________________________________________________
c2690925 2171TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
c04c80e6 2172 //
2173 // Create efficiency grid and calculate efficiency
2174 // of step to step0
2175 //
2176 TString name("eff");
2177 name += step;
2178 name+= step0;
2179 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2180 eff->CalculateEfficiency(step,step0);
2181 return eff;
2182}
c2690925 2183
2184//____________________________________________________________________________
8c1c76e9 2185THnSparse* AliHFEspectrum::GetCharmWeights(){
c2690925 2186
2187 //
2188 // Measured D->e based weighting factors
2189 //
2190
2191 const Int_t nDim=1;
2192 Int_t nBin[nDim] = {44};
2193 const Double_t kPtbound[2] = {0.1, 20.};
2194
2195 Double_t* binEdges[nDim];
2196 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2197
2198 fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
2199 for(Int_t idim = 0; idim < nDim; idim++)
2200 fWeightCharm->SetBinEdges(idim, binEdges[idim]);
8c1c76e9 2201 Double_t weight[44]={
22021.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225};
c2690925 2203 //points
8c1c76e9 2204 Double_t pt[1];
2205 for(int i=0; i<nBin[0]; i++){
2206 pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2207 fWeightCharm->Fill(pt,weight[i]);
c2690925 2208 }
2209 Int_t* ibins[nDim];
2210 for(Int_t ivar = 0; ivar < nDim; ivar++)
2211 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2212 fWeightCharm->SetBinError(ibins[0],0);
2213
2214 return fWeightCharm;
2215}
8c1c76e9 2216
2217//____________________________________________________________________________
2218THnSparse* AliHFEspectrum::GetBeautyIPEff(){
2219 //
2220 // Return beauty electron IP cut efficiency
2221 //
2222
2223 const Int_t nDim=1;
2224 Int_t nBin[nDim] = {44};
2225 const Double_t kPtbound[2] = {0.1, 20.};
2226
2227 Double_t* binEdges[nDim];
2228 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2229
2230 THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin);
2231 for(Int_t idim = 0; idim < nDim; idim++)
2232 ipcut->SetBinEdges(idim, binEdges[idim]);
2233 Double_t pt[1];
2234 Double_t weight;
2235 for(int i=0; i<nBin[0]; i++){
2236 pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2237 weight=0.612*(0.385+0.111*log(pt[0])-0.00435*log(pt[0])*log(pt[0]))*tanh(5.42*pt[0]-1.22); // for 3 sigma cut
2238 //weight=0.786*(0.493+0.0283*log(pt[0])-0.0190*log(pt[0])*log(pt[0]))*tanh(1.84*pt[0]-0.00368); // for 2 sigma cut
2239 ipcut->Fill(pt,weight);
2240 }
2241 Int_t* ibins[nDim];
2242 for(Int_t ivar = 0; ivar < nDim; ivar++)
2243 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2244 ipcut->SetBinError(ibins[0],0);
2245
2246 return ipcut;
2247}
2248
2249//____________________________________________________________________________
2250THnSparse* AliHFEspectrum::GetCharmEff(){
2251 //
2252 // Return charm electron IP cut efficiency
2253 //
2254
2255 const Int_t nDim=1;
2256 Int_t nBin[nDim] = {44};
2257 const Double_t kPtbound[2] = {0.1, 20.};
2258
2259 Double_t* binEdges[nDim];
2260 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2261
2262 THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin);
2263 for(Int_t idim = 0; idim < nDim; idim++)
2264 ipcut->SetBinEdges(idim, binEdges[idim]);
2265 Double_t pt[1];
2266 Double_t weight;
2267 for(int i=0; i<nBin[0]; i++){
2268 pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2269 weight=0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
2270 //weight=0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
2271 ipcut->Fill(pt,weight);
2272 }
2273 Int_t* ibins[nDim];
2274 for(Int_t ivar = 0; ivar < nDim; ivar++)
2275 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2276 ipcut->SetBinError(ibins[0],0);
2277
2278 return ipcut;
2279}
2280
2281//____________________________________________________________________________
2282THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
2283 //
2284 // Return PID x IP cut efficiency
2285 //
2286
2287 const Int_t nDim=1;
2288 Int_t nBin[nDim] = {44};
2289 const Double_t kPtbound[2] = {0.1, 20.};
2290
2291 Double_t* binEdges[nDim];
2292 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2293
2294 THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
2295 for(Int_t idim = 0; idim < nDim; idim++)
2296 pideff->SetBinEdges(idim, binEdges[idim]);
2297
2298 Double_t pt[1];
2299 Double_t weight = 1.0;
2300 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
2301 for(int i=0; i<nBin[0]; i++){
2302 pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2303 Double_t tofeff = 0.9885*(0.8551+0.0070*log(pt[0])-0.0014*log(pt[0])*log(pt[0]))*tanh(0.8884*pt[0]+0.6061); // parameterized TOF PID eff
2304
2305 if(source==0) weight = tofeff*trdtpcPidEfficiency*0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
2306 //if(source==0) weight = tofeff*trdtpcPidEfficiency*0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
2307 //if(source==2) weight = tofeff*trdtpcPidEfficiency*0.5575*(0.3594-0.3051*log(pt[0])+0.1597*log(pt[0])*log(pt[0]))*tanh(8.2436*pt[0]-0.8125);
2308 //if(source==3) weight = tofeff*trdtpcPidEfficiency*0.0937*(0.4449-0.3769*log(pt[0])+0.5031*log(pt[0])*log(pt[0]))*tanh(6.4063*pt[0]+3.1425); // multiply ip cut eff for Dalitz/dielectrons
2309 if(source==2) weight = tofeff*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // multiply ip cut eff for conversion electrons
2310 if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons
2311 printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1));
2312
2313 pideff->Fill(pt,weight);
2314 }
2315 Int_t* ibins[nDim];
2316 for(Int_t ivar = 0; ivar < nDim; ivar++)
2317 ibins[ivar] = new Int_t[nBin[ivar] + 1];
2318 pideff->SetBinError(ibins[0],0);
2319
2320 return pideff;
2321}