2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 // Class for spectrum correction
18 // Subtraction of hadronic background, Unfolding of the data and
19 // Renormalization done here
20 // The following containers have to be set:
21 // - Correction framework container for real data
22 // - Correction framework container for MC (Efficiency Map)
23 // - Correction framework container for background coming from data
24 // - Correction framework container for background coming from MC
27 // Raphaelle Bailhache <R.Bailhache@gsi.de>
28 // Markus Fasel <M.Fasel@gsi.de>
34 #include <TObjArray.h>
41 #include <TGraphErrors.h>
48 #include "AliCFContainer.h"
49 #include "AliCFDataGrid.h"
50 #include "AliCFEffGrid.h"
51 #include "AliCFGridSparse.h"
52 #include "AliCFUnfolding.h"
55 #include "AliHFEspectrum.h"
56 #include "AliHFEcuts.h"
57 #include "AliHFEcontainer.h"
58 #include "AliHFEtools.h"
60 ClassImp(AliHFEspectrum)
62 //____________________________________________________________
63 AliHFEspectrum::AliHFEspectrum(const char *name):
65 fCFContainers(new TObjArray(kDataContainerV0+1)),
66 fTemporaryObjects(NULL),
69 fEfficiencyFunction(NULL),
71 fInclusiveSpectrum(kTRUE),
74 fUnSetCorrelatedErrors(kTRUE),
75 fSetSmoothing(kFALSE),
76 fIPanaHadronBgSubtract(kFALSE),
77 fIPanaCharmBgSubtract(kFALSE),
78 fIPanaConversionBgSubtract(kFALSE),
79 fIPanaNonHFEBgSubtract(kFALSE),
80 fIPParameterizedEff(kFALSE),
82 fBeauty2ndMethod(kFALSE),
83 fIPEffCombinedSamples(kTRUE),
89 fStepBeforeCutsV0(-1),
91 fStepGuessedUnfolding(-1),
92 fNumberOfIterations(5),
93 fChargeChoosen(kAllCharge),
94 fNCentralityBinAtTheEnd(0),
95 fTestCentralityLow(-1),
96 fTestCentralityHigh(-1),
97 fFillMoreCorrelationMatrix(kFALSE),
98 fHadronEffbyIPcut(NULL),
99 fConversionEffbgc(NULL),
101 fBSpectrum2ndMethod(NULL),
102 fkBeauty2ndMethodfilename(""),
106 fWriteToFile(kFALSE),
110 // Default constructor
113 for(Int_t k = 0; k < 20; k++){
116 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
117 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
120 fEfficiencyTOFPIDD[k] = 0;
121 fEfficiencyesdTOFPIDD[k] = 0;
122 fEfficiencyIPCharmD[k] = 0;
123 fEfficiencyIPBeautyD[k] = 0;
124 fEfficiencyIPBeautyesdD[k] = 0;
125 fEfficiencyIPConversionD[k] = 0;
126 fEfficiencyIPNonhfeD[k] = 0;
128 fConversionEff[k] = 0;
132 fEfficiencyCharmSigD[k] = 0;
133 fEfficiencyBeautySigD[k] = 0;
134 fEfficiencyBeautySigesdD[k] = 0;
137 memset(fEtaRange, 0, sizeof(Double_t) * 2);
138 memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
139 memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
140 memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
143 //____________________________________________________________
144 AliHFEspectrum::~AliHFEspectrum(){
148 if(fCFContainers) delete fCFContainers;
149 if(fTemporaryObjects){
150 fTemporaryObjects->Clear();
151 delete fTemporaryObjects;
154 //____________________________________________________________
155 Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
157 // Init what we need for the correction:
159 // Raw spectrum, hadron contamination
160 // MC efficiency maps, correlation matrix
161 // V0 efficiency if wanted
163 // This for a given dimension.
164 // If no inclusive spectrum, then take only efficiency map for beauty electron
165 // and the appropriate correlation matrix
169 if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
174 if(fBeamType==0) kNdim=3;
183 // Get the requested format
186 switch(fNbDimensions){
189 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
191 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
194 AliError("Container with this number of dimensions not foreseen (yet)");
200 // PbPb analysis; centrality as first dimension
201 Int_t nbDimensions = fNbDimensions;
202 fNbDimensions = fNbDimensions + 1;
203 switch(nbDimensions){
208 for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
211 for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
214 AliError("Container with this number of dimensions not foreseen (yet)");
219 // Data container: raw spectrum + hadron contamination
220 AliCFContainer *datacontainer = 0x0;
221 if(fInclusiveSpectrum) {
222 datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
226 datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
228 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
229 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
231 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
232 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
233 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
235 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
236 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
238 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
239 AliCFContainer *mccontaineresd = 0x0;
240 AliCFContainer *mccontaineresdbg = 0x0;
241 AliCFContainer *mccontainermc = 0x0;
242 AliCFContainer *mccontainermcbg = 0x0;
243 AliCFContainer *nonHFEweightedContainer = 0x0;
244 AliCFContainer *convweightedContainer = 0x0;
245 AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
246 AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
248 if(fInclusiveSpectrum) {
249 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
250 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
253 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
254 mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
255 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
256 mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
259 const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
260 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
261 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
262 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
263 nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
264 convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
265 for(Int_t icentr=0;icentr<kNcentr;icentr++)
269 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
270 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
275 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
276 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
278 // if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
280 if(fBeamType == 1)break;
285 nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
286 convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
287 if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
290 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
293 if(!fInclusiveSpectrum) source = 1; //beauty
294 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
295 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
296 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
297 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
298 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
300 // set charm, nonHFE container to estimate BG
301 if(!fInclusiveSpectrum){
303 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
304 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
307 AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
308 SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
309 AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
310 SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
313 SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
316 // MC container: correlation matrix
317 THnSparseF *mccorrelation = 0x0;
318 if(fInclusiveSpectrum) {
319 if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
320 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
321 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
322 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
323 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
325 if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
327 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
328 //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
329 if(!mccorrelation) return kFALSE;
330 Int_t testCentralityLow = fTestCentralityLow;
331 Int_t testCentralityHigh = fTestCentralityHigh;
332 if(fFillMoreCorrelationMatrix) {
333 testCentralityLow = fTestCentralityLow-1;
334 testCentralityHigh = fTestCentralityHigh+1;
336 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
337 if(!mccorrelationD) {
338 printf("No correlation\n");
341 SetCorrelation(mccorrelationD);
343 // V0 container Electron, pt eta phi
345 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
346 if(!containerV0) return kFALSE;
347 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
348 if(!containerV0Electron) return kFALSE;
349 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
355 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
356 contaminationspectrum->SetName("contaminationspectrum");
357 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
358 ccontaminationspectrum->Divide(3,1);
359 ccontaminationspectrum->cd(1);
360 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
361 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
362 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
363 contaminationspectrum2dpteta->SetStats(0);
364 contaminationspectrum2dpteta->SetTitle("");
365 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
366 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
367 contaminationspectrum2dptphi->SetStats(0);
368 contaminationspectrum2dptphi->SetTitle("");
369 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
370 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
371 contaminationspectrum2detaphi->SetStats(0);
372 contaminationspectrum2detaphi->SetTitle("");
373 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
374 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
375 contaminationspectrum2dptphi->Draw("colz");
376 ccontaminationspectrum->cd(2);
377 contaminationspectrum2dpteta->Draw("colz");
378 ccontaminationspectrum->cd(3);
379 contaminationspectrum2detaphi->Draw("colz");
381 TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
385 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
386 ptcorrelation->Draw("colz");
389 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
390 ptcorrelation->Draw("colz");
394 ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
395 ccorrelation->SaveAs("correlationMatrix.eps");
399 TFile *file = TFile::Open("tests.root","recreate");
400 datacontainerD->Write("data");
401 mccontainermcD->Write("mcefficiency");
402 mccorrelationD->Write("correlationmatrix");
409 //____________________________________________________________
410 void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
412 // get spectrum for beauty 2nd method
415 TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
416 fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
420 //____________________________________________________________
421 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
423 // Correct the spectrum for efficiency and unfolding
424 // with both method and compare
427 gStyle->SetPalette(1);
428 gStyle->SetOptStat(1111);
429 gStyle->SetPadBorderMode(0);
430 gStyle->SetCanvasColor(10);
431 gStyle->SetPadLeftMargin(0.13);
432 gStyle->SetPadRightMargin(0.13);
434 printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
436 ///////////////////////////
437 // Check initialization
438 ///////////////////////////
440 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
441 AliInfo("You have to init before");
445 if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
446 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
450 SetNumberOfIteration(10);
451 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
453 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
454 //////////////////////////////////
455 // Subtract hadron background
456 /////////////////////////////////
457 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
458 if(subtractcontamination) {
459 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
460 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
463 ////////////////////////////////////////////////
464 // Correct for TPC efficiency from V0
465 ///////////////////////////////////////////////
466 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
467 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
468 if(dataContainerV0){printf("Got the V0 container\n");
469 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
470 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
473 //////////////////////////////////////////////////////////////////////////////
474 // Correct for efficiency parametrized (if TPC efficiency is parametrized)
475 /////////////////////////////////////////////////////////////////////////////
476 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
477 if(fEfficiencyFunction){
478 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
479 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
485 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
487 printf("Unfolded failed\n");
490 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
491 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
492 if(!correctedspectrum){
493 AliError("No corrected spectrum\n");
496 if(!residualspectrum){
497 AliError("No residul spectrum\n");
501 /////////////////////
504 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
510 if(fDebugLevel > 0.0) {
513 if(fBeamType==0) ptpr=0;
514 if(fBeamType==1) ptpr=1;
516 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
517 ccorrected->Divide(2,1);
520 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
521 correctedspectrumD->SetTitle("");
522 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
523 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
524 correctedspectrumD->SetMarkerStyle(26);
525 correctedspectrumD->SetMarkerColor(kBlue);
526 correctedspectrumD->SetLineColor(kBlue);
527 correctedspectrumD->Draw("AP");
528 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
529 alltogetherspectrumD->SetTitle("");
530 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
531 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
532 alltogetherspectrumD->SetMarkerStyle(25);
533 alltogetherspectrumD->SetMarkerColor(kBlack);
534 alltogetherspectrumD->SetLineColor(kBlack);
535 alltogetherspectrumD->Draw("P");
536 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
537 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
538 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
539 legcorrected->Draw("same");
541 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
542 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
543 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
544 ratiocorrected->SetName("ratiocorrected");
545 ratiocorrected->SetTitle("");
546 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
547 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
548 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
549 ratiocorrected->SetStats(0);
550 ratiocorrected->Draw();
551 if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
553 //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
554 //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
555 //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
556 //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
558 TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
559 TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
560 TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
561 TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
565 TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
566 ccorrectedallspectra->Divide(2,1);
567 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
568 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
570 THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
571 TAxis *cenaxisa = sparsesu->GetAxis(0);
572 THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
573 TAxis *cenaxisb = sparsed->GetAxis(0);
574 Int_t nbbin = cenaxisb->GetNbins();
575 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
576 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
577 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
578 TString titlee("corrected_centrality_bin_");
580 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
582 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
584 TString titleec("corrected_check_projection_bin_");
586 titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
588 titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
590 TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
591 titleenameunotnormalized += "[";
592 titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
593 titleenameunotnormalized += "_";
594 titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
595 titleenameunotnormalized += "[";
596 TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
597 titleenameunormalized += "[";
598 titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
599 titleenameunormalized += "_";
600 titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
601 titleenameunormalized += "[";
602 TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
603 titleenamednotnormalized += "[";
604 titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
605 titleenamednotnormalized += "_";
606 titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
607 titleenamednotnormalized += "[";
608 TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
609 titleenamednormalized += "[";
610 titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
611 titleenamednormalized += "_";
612 titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
613 titleenamednormalized += "[";
615 for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
616 printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
617 nbEvents += fNEvents[k];
619 Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
620 Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
621 printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
622 Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
623 Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
624 printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
625 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
626 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
627 TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
628 ccorrectedcheck->cd(1);
629 TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
630 TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
632 aftersdc->Draw("same");
633 TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
634 ccorrectede->Divide(2,1);
637 TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
638 CorrectFromTheWidth(aftersu);
639 aftersu->SetName((const char*)titleenameunotnormalized);
640 unfoldingspectrac[binc] = *aftersu;
642 TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
643 aftersun->SetTitle("");
644 aftersun->GetYaxis()->SetTitleOffset(1.5);
645 aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
646 aftersun->SetMarkerStyle(26);
647 aftersun->SetMarkerColor(kBlue);
648 aftersun->SetLineColor(kBlue);
649 aftersun->Draw("AP");
650 aftersun->SetName((const char*)titleenameunormalized);
651 unfoldingspectracn[binc] = *aftersun;
653 TH1D *aftersd = (TH1D *) sparsed->Projection(1);
654 CorrectFromTheWidth(aftersd);
655 aftersd->SetName((const char*)titleenamednotnormalized);
656 correctedspectrac[binc] = *aftersd;
658 TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
659 aftersdn->SetTitle("");
660 aftersdn->GetYaxis()->SetTitleOffset(1.5);
661 aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
662 aftersdn->SetMarkerStyle(25);
663 aftersdn->SetMarkerColor(kBlack);
664 aftersdn->SetLineColor(kBlack);
666 aftersdn->SetName((const char*)titleenamednormalized);
667 correctedspectracn[binc] = *aftersdn;
668 TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
669 legcorrectedud->AddEntry(aftersun,"Corrected","p");
670 legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
671 legcorrectedud->Draw("same");
672 ccorrectedallspectra->cd(1);
674 TH1D *aftersunn = (TH1D *) aftersun->Clone();
675 aftersunn->SetMarkerStyle(stylee[binc]);
676 aftersunn->SetMarkerColor(colorr[binc]);
677 if(binc==0) aftersunn->Draw("AP");
678 else aftersunn->Draw("P");
679 legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
680 ccorrectedallspectra->cd(2);
682 TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
683 aftersdnn->SetMarkerStyle(stylee[binc]);
684 aftersdnn->SetMarkerColor(colorr[binc]);
685 if(binc==0) aftersdnn->Draw("AP");
686 else aftersdnn->Draw("P");
687 legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
689 TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
690 TString titleee("ratiocorrected_bin_");
692 ratiocorrectedbinc->SetName((const char*) titleee);
693 ratiocorrectedbinc->SetTitle("");
694 ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
695 ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
696 ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
697 ratiocorrectedbinc->SetStats(0);
698 ratiocorrectedbinc->Draw();
701 ccorrectedallspectra->cd(1);
702 legtotal->Draw("same");
703 ccorrectedallspectra->cd(2);
704 legtotalg->Draw("same");
706 cenaxisa->SetRange(0,nbbin);
707 cenaxisb->SetRange(0,nbbin);
708 if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
711 // Dump to file if needed
713 TFile *out = new TFile("finalSpectrum.root","recreate");
714 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
715 correctedspectrumD->Write();
716 alltogetherspectrumD->SetName("AlltogetherSpectrum");
717 alltogetherspectrumD->Write();
718 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
719 ratiocorrected->Write();
720 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
721 correctedspectrum->Write();
722 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
723 alltogetherCorrection->Write();
724 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
725 unfoldingspectrac[binc].Write();
726 unfoldingspectracn[binc].Write();
727 correctedspectrac[binc].Write();
728 correctedspectracn[binc].Write();
730 out->Close(); delete out;
733 if (unfoldingspectrac) delete[] unfoldingspectrac;
734 if (unfoldingspectracn) delete[] unfoldingspectracn;
735 if (correctedspectrac) delete[] correctedspectrac;
736 if (correctedspectracn) delete[] correctedspectracn;
743 //____________________________________________________________
744 Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
746 // Correct the spectrum for efficiency and unfolding for beauty analysis
747 // with both method and compare
750 gStyle->SetPalette(1);
751 gStyle->SetOptStat(1111);
752 gStyle->SetPadBorderMode(0);
753 gStyle->SetCanvasColor(10);
754 gStyle->SetPadLeftMargin(0.13);
755 gStyle->SetPadRightMargin(0.13);
757 ///////////////////////////
758 // Check initialization
759 ///////////////////////////
761 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
762 AliInfo("You have to init before");
766 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
767 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
771 SetNumberOfIteration(10);
772 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
774 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
775 //////////////////////////////////
776 // Subtract hadron background
777 /////////////////////////////////
778 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
779 AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
780 TGraphErrors *gNormalizedRawSpectrum = 0x0;
781 if(subtractcontamination) {
782 if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
783 else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
784 unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
785 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
786 gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
789 printf("after normalize getting IP \n");
791 /////////////////////////////////////////////////////////////////////////////////////////
792 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
793 /////////////////////////////////////////////////////////////////////////////////////////
795 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
796 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
797 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
799 if(fEfficiencyFunction){
800 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
801 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
803 else if(dataContainerV0){
804 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
805 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
813 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
815 printf("Unfolded failed\n");
818 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
819 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
820 if(!correctedspectrum){
821 AliError("No corrected spectrum\n");
824 if(!residualspectrum){
825 AliError("No residual spectrum\n");
829 /////////////////////
833 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
840 if(fDebugLevel > 0.0) {
843 if(fBeamType==0) ptpr=0;
844 if(fBeamType==1) ptpr=1;
846 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
847 ccorrected->Divide(2,1);
850 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
851 correctedspectrumD->SetTitle("");
852 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
853 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
854 correctedspectrumD->SetMarkerStyle(26);
855 correctedspectrumD->SetMarkerColor(kBlue);
856 correctedspectrumD->SetLineColor(kBlue);
857 correctedspectrumD->Draw("AP");
858 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
859 alltogetherspectrumD->SetTitle("");
860 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
861 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
862 alltogetherspectrumD->SetMarkerStyle(25);
863 alltogetherspectrumD->SetMarkerColor(kBlack);
864 alltogetherspectrumD->SetLineColor(kBlack);
865 alltogetherspectrumD->Draw("P");
866 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
867 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
868 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
869 legcorrected->Draw("same");
871 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
872 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
873 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
874 ratiocorrected->SetName("ratiocorrected");
875 ratiocorrected->SetTitle("");
876 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
877 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
878 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
879 ratiocorrected->SetStats(0);
880 ratiocorrected->Draw();
881 if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
885 CalculateNonHFEsyst(0);
889 // Dump to file if needed
892 // to do centrality dependent
895 out = new TFile("finalSpectrum.root","recreate");
898 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
899 correctedspectrumD->Write();
900 alltogetherspectrumD->SetName("AlltogetherSpectrum");
901 alltogetherspectrumD->Write();
902 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
903 ratiocorrected->Write();
905 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
906 correctedspectrum->Write();
907 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
908 alltogetherCorrection->Write();
910 if(unnormalizedRawSpectrum) {
911 unnormalizedRawSpectrum->SetName("beautyAfterIP");
912 unnormalizedRawSpectrum->Write();
915 if(gNormalizedRawSpectrum){
916 gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
917 gNormalizedRawSpectrum->Write();
922 fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
923 fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
924 fEfficiencyCharmSigD[countpp]->Write();
925 fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
926 fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
927 fEfficiencyBeautySigD[countpp]->Write();
928 fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
929 fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
930 fCharmEff[countpp]->Write();
931 fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
932 fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
933 fBeautyEff[countpp]->Write();
934 fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
935 fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
936 fConversionEff[countpp]->Write();
937 fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
938 fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
939 fNonHFEEff[countpp]->Write();
944 TGraphErrors* correctedspectrumDc[kCentrality];
945 TGraphErrors* alltogetherspectrumDc[kCentrality];
946 for(Int_t i=0;i<kCentrality-2;i++)
948 correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
949 correctedspectrumDc[i] = Normalize(correctedspectrum,i);
950 if(correctedspectrumDc[i]){
951 correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
952 correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
953 correctedspectrumDc[i]->Write();
955 alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
956 alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
957 if(alltogetherspectrumDc[i]){
958 alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
959 alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
960 alltogetherspectrumDc[i]->Write();
963 TH1D *centrcrosscheck = correctedspectrum->Projection(0);
964 centrcrosscheck->SetTitle(Form("centrality_%i",i));
965 centrcrosscheck->SetName(Form("centrality_%i",i));
966 centrcrosscheck->Write();
968 TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
969 TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
971 TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0);
972 centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
973 centrcrosscheck2->SetName(Form("centrality2_%i",i));
974 centrcrosscheck2->Write();
976 TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
977 ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
978 ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
979 ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
980 ratiocorrectedc->Write();
982 fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
983 fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
984 fEfficiencyCharmSigD[i]->Write();
985 fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
986 fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
987 fEfficiencyBeautySigD[i]->Write();
988 fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
989 fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
990 fCharmEff[i]->Write();
991 fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
992 fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
993 fBeautyEff[i]->Write();
994 fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
995 fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
996 fConversionEff[i]->Write();
997 fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
998 fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
999 fNonHFEEff[i]->Write();
1012 //____________________________________________________________
1013 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
1015 // Apply background subtraction
1032 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1034 AliError("Data Container not available");
1037 printf("Step data: %d\n",fStepData);
1038 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
1040 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
1041 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
1044 // Background Estimate
1045 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
1046 if(!backgroundContainer){
1047 AliError("MC background container not found");
1051 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
1052 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
1054 if(!fInclusiveSpectrum){
1055 //Background subtraction for IP analysis
1057 TH1D *incElecCent[kCentrality-1];
1058 TH1D *charmCent[kCentrality-1];
1059 TH1D *convCent[kCentrality-1];
1060 TH1D *nonHFECent[kCentrality-1];
1061 TH1D *subtractedCent[kCentrality-1];
1062 TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
1063 CorrectFromTheWidth(measuredTH1Draw);
1065 THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1066 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1067 sparseIncElec->GetAxis(0)->SetRange(icent,icent);
1068 incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
1069 CorrectFromTheWidth(incElecCent[icent-1]);
1072 TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
1074 rawspectra->SetLogy();
1075 gStyle->SetOptStat(0);
1076 TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
1077 measuredTH1Draw->SetMarkerStyle(20);
1078 measuredTH1Draw->Draw();
1079 measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
1080 lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
1082 Int_t* bins=new Int_t[2];
1083 if(fIPanaHadronBgSubtract){
1084 // Hadron background
1085 printf("Hadron background for IP analysis subtracted!\n");
1089 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1090 bins[0]=htemp->GetNbinsX();
1094 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1095 bins[0]=htemp->GetNbinsX();
1096 htemp = (TH1D *) fHadronEffbyIPcut->Projection(1);
1097 bins[1]=htemp->GetNbinsX();
1099 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
1100 hbgContainer->SetGrid(fHadronEffbyIPcut);
1101 backgroundGrid->Multiply(hbgContainer,1);
1102 // draw raw hadron bg spectra
1103 TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
1104 CorrectFromTheWidth(hadronbg);
1105 hadronbg->SetMarkerColor(7);
1106 hadronbg->SetMarkerStyle(20);
1108 hadronbg->Draw("samep");
1109 lRaw->AddEntry(hadronbg,"hadrons");
1110 // subtract hadron contamination
1111 spectrumSubtracted->Add(backgroundGrid,-1.0);
1113 if(fIPanaCharmBgSubtract){
1115 printf("Charm background for IP analysis subtracted!\n");
1116 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
1117 // draw charm bg spectra
1118 TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
1119 CorrectFromTheWidth(charmbg);
1120 charmbg->SetMarkerColor(3);
1121 charmbg->SetMarkerStyle(20);
1123 charmbg->Draw("samep");
1124 lRaw->AddEntry(charmbg,"charm elecs");
1125 // subtract charm background
1126 spectrumSubtracted->Add(charmbgContainer,-1.0);
1128 THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
1129 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1130 sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
1131 charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
1132 CorrectFromTheWidth(charmCent[icent-1]);
1136 if(fIPanaConversionBgSubtract){
1137 // Conversion background
1138 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
1139 // draw conversion bg spectra
1140 TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
1141 CorrectFromTheWidth(conversionbg);
1142 conversionbg->SetMarkerColor(4);
1143 conversionbg->SetMarkerStyle(20);
1145 conversionbg->Draw("samep");
1146 lRaw->AddEntry(conversionbg,"conversion elecs");
1147 // subtract conversion background
1148 spectrumSubtracted->Add(conversionbgContainer,-1.0);
1150 THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
1151 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1152 sparseconvElec->GetAxis(0)->SetRange(icent,icent);
1153 convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
1154 CorrectFromTheWidth(convCent[icent-1]);
1158 if(fIPanaNonHFEBgSubtract){
1159 // NonHFE background
1160 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
1161 // draw Dalitz/dielectron bg spectra
1162 TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
1163 CorrectFromTheWidth(nonhfebg);
1164 nonhfebg->SetMarkerColor(6);
1165 nonhfebg->SetMarkerStyle(20);
1167 nonhfebg->Draw("samep");
1168 lRaw->AddEntry(nonhfebg,"non-HF elecs");
1169 // subtract Dalitz/dielectron background
1170 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
1172 THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
1173 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1174 sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
1175 nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
1176 CorrectFromTheWidth(nonHFECent[icent-1]);
1181 TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
1182 CorrectFromTheWidth(rawbgsubtracted);
1183 rawbgsubtracted->SetMarkerStyle(24);
1185 lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
1186 rawbgsubtracted->Draw("samep");
1189 //rawspectra->SaveAs("rawspectra.eps");
1192 THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1193 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1194 sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
1195 subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
1196 CorrectFromTheWidth(subtractedCent[icent-1]);
1199 TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
1200 TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
1201 centRaw->Divide(3,3);
1202 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1206 incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
1207 incElecCent[icent-1]->Draw("p");
1208 incElecCent[icent-1]->SetMarkerColor(1);
1209 incElecCent[icent-1]->SetMarkerStyle(20);
1210 charmCent[icent-1]->Draw("samep");
1211 charmCent[icent-1]->SetMarkerColor(3);
1212 charmCent[icent-1]->SetMarkerStyle(20);
1213 convCent[icent-1]->Draw("samep");
1214 convCent[icent-1]->SetMarkerColor(4);
1215 convCent[icent-1]->SetMarkerStyle(20);
1216 nonHFECent[icent-1]->Draw("samep");
1217 nonHFECent[icent-1]->SetMarkerColor(6);
1218 nonHFECent[icent-1]->SetMarkerStyle(20);
1219 subtractedCent[icent-1]->Draw("samep");
1220 subtractedCent[icent-1]->SetMarkerStyle(24);
1222 lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
1223 lCentRaw->AddEntry(charmCent[0],"charm elecs");
1224 lCentRaw->AddEntry(convCent[0],"conversion elecs");
1225 lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
1226 lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
1227 lCentRaw->Draw("SAME");
1237 spectrumSubtracted->Add(backgroundGrid,-1.0);
1241 if(fBackground) delete fBackground;
1242 fBackground = backgroundGrid;
1243 } else delete backgroundGrid;
1246 if(fDebugLevel > 0) {
1249 if(fBeamType==0) ptprd=0;
1250 if(fBeamType==1) ptprd=1;
1252 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
1253 cbackgroundsubtraction->Divide(3,1);
1254 cbackgroundsubtraction->cd(1);
1256 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
1257 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
1258 CorrectFromTheWidth(measuredTH1Daftersubstraction);
1259 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
1260 measuredTH1Daftersubstraction->SetStats(0);
1261 measuredTH1Daftersubstraction->SetTitle("");
1262 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1263 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1264 measuredTH1Daftersubstraction->SetMarkerStyle(25);
1265 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
1266 measuredTH1Daftersubstraction->SetLineColor(kBlack);
1267 measuredTH1Dbeforesubstraction->SetStats(0);
1268 measuredTH1Dbeforesubstraction->SetTitle("");
1269 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1270 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1271 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
1272 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
1273 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
1274 measuredTH1Daftersubstraction->Draw();
1275 measuredTH1Dbeforesubstraction->Draw("same");
1276 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
1277 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
1278 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
1279 legsubstraction->Draw("same");
1280 cbackgroundsubtraction->cd(2);
1282 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1283 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1284 ratiomeasuredcontamination->SetTitle("");
1285 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1286 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1287 ratiomeasuredcontamination->Sumw2();
1288 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1289 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1290 ratiomeasuredcontamination->SetStats(0);
1291 ratiomeasuredcontamination->SetMarkerStyle(26);
1292 ratiomeasuredcontamination->SetMarkerColor(kBlack);
1293 ratiomeasuredcontamination->SetLineColor(kBlack);
1294 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1295 ratiomeasuredcontamination->SetBinError(k+1,0.0);
1297 ratiomeasuredcontamination->Draw("P");
1298 cbackgroundsubtraction->cd(3);
1299 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
1300 CorrectFromTheWidth(measuredTH1background);
1301 measuredTH1background->SetStats(0);
1302 measuredTH1background->SetTitle("");
1303 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1304 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1305 measuredTH1background->SetMarkerStyle(26);
1306 measuredTH1background->SetMarkerColor(kBlack);
1307 measuredTH1background->SetLineColor(kBlack);
1308 measuredTH1background->Draw();
1309 if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
1313 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1314 cbackgrounde->Divide(2,1);
1315 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1316 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1318 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1319 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1320 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1321 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1322 Int_t nbbin = cenaxisb->GetNbins();
1323 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1324 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1325 for(Int_t binc = 0; binc < nbbin; binc++){
1326 TString titlee("BackgroundSubtraction_centrality_bin_");
1328 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1329 cbackground->Divide(2,1);
1332 cenaxisa->SetRange(binc+1,binc+1);
1333 cenaxisb->SetRange(binc+1,binc+1);
1334 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1335 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1336 CorrectFromTheWidth(aftersubstraction);
1337 CorrectFromTheWidth(beforesubstraction);
1338 aftersubstraction->SetStats(0);
1339 aftersubstraction->SetTitle((const char*)titlee);
1340 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1341 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1342 aftersubstraction->SetMarkerStyle(25);
1343 aftersubstraction->SetMarkerColor(kBlack);
1344 aftersubstraction->SetLineColor(kBlack);
1345 beforesubstraction->SetStats(0);
1346 beforesubstraction->SetTitle((const char*)titlee);
1347 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1348 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1349 beforesubstraction->SetMarkerStyle(24);
1350 beforesubstraction->SetMarkerColor(kBlue);
1351 beforesubstraction->SetLineColor(kBlue);
1352 aftersubstraction->Draw();
1353 beforesubstraction->Draw("same");
1354 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1355 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1356 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1358 cbackgrounde->cd(1);
1360 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1361 aftersubtractionn->SetMarkerStyle(stylee[binc]);
1362 aftersubtractionn->SetMarkerColor(colorr[binc]);
1363 if(binc==0) aftersubtractionn->Draw();
1364 else aftersubtractionn->Draw("same");
1365 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1366 cbackgrounde->cd(2);
1368 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1369 aftersubtractionng->SetMarkerStyle(stylee[binc]);
1370 aftersubtractionng->SetMarkerColor(colorr[binc]);
1371 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1372 if(binc==0) aftersubtractionng->Draw();
1373 else aftersubtractionng->Draw("same");
1374 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1376 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1377 ratiocontamination->SetName("ratiocontamination");
1378 ratiocontamination->SetTitle((const char*)titlee);
1379 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1380 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1381 ratiocontamination->Add(aftersubstraction,-1.0);
1382 ratiocontamination->Divide(beforesubstraction);
1383 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1384 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1385 ratiocontamination->SetBinError(nbinpt+1,0.0);
1387 ratiocontamination->SetStats(0);
1388 ratiocontamination->SetMarkerStyle(26);
1389 ratiocontamination->SetMarkerColor(kBlack);
1390 ratiocontamination->SetLineColor(kBlack);
1391 ratiocontamination->Draw("P");
1394 cbackgrounde->cd(1);
1395 legtotal->Draw("same");
1396 cbackgrounde->cd(2);
1397 legtotalg->Draw("same");
1399 cenaxisa->SetRange(0,nbbin);
1400 cenaxisb->SetRange(0,nbbin);
1401 if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
1405 return spectrumSubtracted;
1408 //____________________________________________________________
1409 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1411 // calculate charm background
1426 if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
1428 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1430 AliError("MC Container not available");
1435 AliError("No Correlation map available");
1439 AliCFDataGrid *charmBackgroundGrid= 0x0;
1440 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1442 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1443 Int_t* bins=new Int_t[2];
1444 bins[0]=charmbgaftertofpid->GetNbinsX();
1449 charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
1450 bins[1]=charmbgaftertofpid->GetNbinsX();
1454 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
1455 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1456 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
1458 charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
1460 Int_t *nBinpp=new Int_t[1];
1461 Int_t* binspp=new Int_t[1];
1462 binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
1464 Int_t *nBinPbPb=new Int_t[2];
1465 Int_t* binsPbPb=new Int_t[2];
1466 binsPbPb[1]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
1469 Int_t looppt=binspp[0];
1470 if(fBeamType==1) looppt=binsPbPb[1];
1472 for(Long_t iBin=1; iBin<= looppt;iBin++){
1476 charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
1477 charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
1481 // loop over centrality
1482 for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1485 Double_t evtnormPbPb=0;
1486 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1487 charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1488 charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1493 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
1495 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
1496 weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1497 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
1498 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1499 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
1501 // Efficiency (set efficiency to 1 for only folding)
1502 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1503 efficiencyD->CalculateEfficiency(0,0);
1506 if(fBeamType==0)nDim = 1;
1507 if(fBeamType==1)nDim = 2;
1508 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1509 folding.SetMaxNumberOfIterations(1);
1513 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1514 THnSparse* result=(THnSparse*)result1->Clone();
1515 charmBackgroundGrid->SetGrid(result);
1516 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
1518 //Charm background evaluation plots
1520 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1521 cCharmBgEval->Divide(3,1);
1523 cCharmBgEval->cd(1);
1525 if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
1528 Double_t evtnormPbPb=0;
1529 for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1531 if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1533 charmbgaftertofpid->Scale(evtnormPbPb);
1536 CorrectFromTheWidth(charmbgaftertofpid);
1537 charmbgaftertofpid->SetMarkerStyle(25);
1538 charmbgaftertofpid->Draw("p");
1539 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1540 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1543 CorrectFromTheWidth(charmbgafteripcut);
1544 charmbgafteripcut->SetMarkerStyle(24);
1545 charmbgafteripcut->Draw("samep");
1547 CorrectFromTheWidth(charmbgafterweight);
1548 charmbgafterweight->SetMarkerStyle(24);
1549 charmbgafterweight->SetMarkerColor(4);
1550 charmbgafterweight->Draw("samep");
1552 CorrectFromTheWidth(charmbgafterfolding);
1553 charmbgafterfolding->SetMarkerStyle(24);
1554 charmbgafterfolding->SetMarkerColor(2);
1555 charmbgafterfolding->Draw("samep");
1557 cCharmBgEval->cd(2);
1558 parametrizedcharmpidipeff->SetMarkerStyle(24);
1559 parametrizedcharmpidipeff->Draw("p");
1560 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1562 cCharmBgEval->cd(3);
1563 charmweightingfc->SetMarkerStyle(24);
1564 charmweightingfc->Draw("p");
1565 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1566 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1568 cCharmBgEval->cd(1);
1569 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1570 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1571 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1572 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1573 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1574 legcharmbg->Draw("same");
1576 cCharmBgEval->cd(2);
1577 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1578 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1579 legcharmbg2->Draw("same");
1581 CorrectStatErr(charmBackgroundGrid);
1582 if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
1590 if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
1592 return charmBackgroundGrid;
1595 //____________________________________________________________
1596 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1598 // calculate conversion background
1601 Double_t evtnorm[1] = {0.0};
1602 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1603 printf("check event!!! %lf \n",evtnorm[0]);
1605 AliCFContainer *backgroundContainer = 0x0;
1608 backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
1609 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1610 backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
1614 // Background Estimate
1615 backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1617 if(!backgroundContainer){
1618 AliError("MC background container not found");
1622 Int_t stepbackground = 3;
1624 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1625 Int_t *nBinpp=new Int_t[1];
1626 Int_t* binspp=new Int_t[1];
1627 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1629 Int_t *nBinPbPb=new Int_t[2];
1630 Int_t* binsPbPb=new Int_t[2];
1631 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1634 Int_t looppt=binspp[0];
1635 if(fBeamType==1) looppt=binsPbPb[1];
1637 for(Long_t iBin=1; iBin<= looppt;iBin++){
1641 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1642 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1646 // loop over centrality
1647 for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1650 Double_t evtnormPbPb=0;
1651 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1652 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1653 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1657 //end of workaround for statistical errors
1659 AliCFDataGrid *weightedConversionContainer;
1660 if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
1661 else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
1662 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1663 backgroundGrid->Multiply(weightedConversionContainer,1.0);
1670 return backgroundGrid;
1674 //____________________________________________________________
1675 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1677 // calculate non-HFE background
1680 Double_t evtnorm[1] = {0.0};
1681 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1682 printf("check event!!! %lf \n",evtnorm[0]);
1684 AliCFContainer *backgroundContainer = 0x0;
1686 backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
1687 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1689 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
1691 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
1695 // Background Estimate
1696 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1698 if(!backgroundContainer){
1699 AliError("MC background container not found");
1704 Int_t stepbackground = 3;
1706 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1707 Int_t *nBinpp=new Int_t[1];
1708 Int_t* binspp=new Int_t[1];
1709 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1711 Int_t *nBinPbPb=new Int_t[2];
1712 Int_t* binsPbPb=new Int_t[2];
1713 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1716 Int_t looppt=binspp[0];
1717 if(fBeamType==1) looppt=binsPbPb[1];
1720 for(Long_t iBin=1; iBin<= looppt;iBin++){
1724 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1725 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1729 for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1732 Double_t evtnormPbPb=0;
1733 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1734 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1735 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1739 //end of workaround for statistical errors
1740 AliCFDataGrid *weightedNonHFEContainer;
1741 if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1742 else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
1743 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1744 backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1751 return backgroundGrid;
1754 //____________________________________________________________
1755 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1758 // Apply TPC pid efficiency correction from parametrisation
1761 // Data in the right format
1762 AliCFDataGrid *dataGrid = 0x0;
1764 dataGrid = bgsubpectrum;
1768 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1770 AliError("Data Container not available");
1773 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1775 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1776 result->SetName("ParametrizedEfficiencyBefore");
1777 THnSparse *h = result->GetGrid();
1778 Int_t nbdimensions = h->GetNdimensions();
1779 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1780 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1782 AliError("Data Container not available");
1785 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1786 dataContainerbis->Add(dataContainerbis,-1.0);
1789 Int_t* coord = new Int_t[nbdimensions];
1790 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1791 Double_t* points = new Double_t[nbdimensions];
1793 ULong64_t nEntries = h->GetNbins();
1794 for (ULong64_t i = 0; i < nEntries; ++i) {
1796 Double_t value = h->GetBinContent(i, coord);
1797 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1798 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1800 // Get the bin co-ordinates given an coord
1801 for (Int_t j = 0; j < nbdimensions; ++j)
1802 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1804 if (!fEfficiencyFunction->IsInside(points))
1806 TF1::RejectPoint(kFALSE);
1808 // Evaulate function at points
1809 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1810 //printf("Value efficiency is %f\n",valueEfficiency);
1812 if(valueEfficiency > 0.0) {
1813 h->SetBinContent(coord,value/valueEfficiency);
1814 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1816 Double_t error = h->GetBinError(i);
1817 h->SetBinError(coord,error/valueEfficiency);
1818 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1826 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1828 if(fDebugLevel > 0) {
1830 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1831 cEfficiencyParametrized->Divide(2,1);
1832 cEfficiencyParametrized->cd(1);
1833 TH1D *afterE = (TH1D *) resultt->Project(0);
1834 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1835 CorrectFromTheWidth(afterE);
1836 CorrectFromTheWidth(beforeE);
1837 afterE->SetStats(0);
1838 afterE->SetTitle("");
1839 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1840 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1841 afterE->SetMarkerStyle(25);
1842 afterE->SetMarkerColor(kBlack);
1843 afterE->SetLineColor(kBlack);
1844 beforeE->SetStats(0);
1845 beforeE->SetTitle("");
1846 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1847 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1848 beforeE->SetMarkerStyle(24);
1849 beforeE->SetMarkerColor(kBlue);
1850 beforeE->SetLineColor(kBlue);
1853 beforeE->Draw("same");
1854 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1855 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1856 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1857 legefficiencyparametrized->Draw("same");
1858 cEfficiencyParametrized->cd(2);
1859 fEfficiencyFunction->Draw();
1860 //cEfficiencyParametrized->cd(3);
1861 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1862 //ratioefficiency->Divide(afterE);
1863 //ratioefficiency->Draw();
1865 if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
1872 //____________________________________________________________
1873 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1876 // Apply TPC pid efficiency correction from V0
1879 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1881 AliError("V0 Container not available");
1886 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1887 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1889 // Data in the right format
1890 AliCFDataGrid *dataGrid = 0x0;
1892 dataGrid = bgsubpectrum;
1896 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1898 AliError("Data Container not available");
1902 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1906 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1907 result->ApplyEffCorrection(*efficiencyD);
1909 if(fDebugLevel > 0) {
1912 if(fBeamType==0) ptpr=0;
1913 if(fBeamType==1) ptpr=1;
1915 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1916 cV0Efficiency->Divide(2,1);
1917 cV0Efficiency->cd(1);
1918 TH1D *afterE = (TH1D *) result->Project(ptpr);
1919 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1920 afterE->SetStats(0);
1921 afterE->SetTitle("");
1922 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1923 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1924 afterE->SetMarkerStyle(25);
1925 afterE->SetMarkerColor(kBlack);
1926 afterE->SetLineColor(kBlack);
1927 beforeE->SetStats(0);
1928 beforeE->SetTitle("");
1929 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1930 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1931 beforeE->SetMarkerStyle(24);
1932 beforeE->SetMarkerColor(kBlue);
1933 beforeE->SetLineColor(kBlue);
1935 beforeE->Draw("same");
1936 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1937 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1938 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1939 legV0efficiency->Draw("same");
1940 cV0Efficiency->cd(2);
1941 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1942 efficiencyDproj->SetTitle("");
1943 efficiencyDproj->SetStats(0);
1944 efficiencyDproj->SetMarkerStyle(25);
1945 efficiencyDproj->Draw();
1949 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1950 cV0Efficiencye->Divide(2,1);
1951 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1952 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1954 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1955 TAxis *cenaxisa = sparseafter->GetAxis(0);
1956 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1957 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1958 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1959 TAxis *cenaxisc = efficiencya->GetAxis(0);
1960 Int_t nbbin = cenaxisb->GetNbins();
1961 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1962 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1963 for(Int_t binc = 0; binc < nbbin; binc++){
1964 TString titlee("V0Efficiency_centrality_bin_");
1966 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1967 ccV0Efficiency->Divide(2,1);
1968 ccV0Efficiency->cd(1);
1970 cenaxisa->SetRange(binc+1,binc+1);
1971 cenaxisb->SetRange(binc+1,binc+1);
1972 cenaxisc->SetRange(binc+1,binc+1);
1973 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1974 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1975 CorrectFromTheWidth(aftere);
1976 CorrectFromTheWidth(beforee);
1977 aftere->SetStats(0);
1978 aftere->SetTitle((const char*)titlee);
1979 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1980 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1981 aftere->SetMarkerStyle(25);
1982 aftere->SetMarkerColor(kBlack);
1983 aftere->SetLineColor(kBlack);
1984 beforee->SetStats(0);
1985 beforee->SetTitle((const char*)titlee);
1986 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1987 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1988 beforee->SetMarkerStyle(24);
1989 beforee->SetMarkerColor(kBlue);
1990 beforee->SetLineColor(kBlue);
1992 beforee->Draw("same");
1993 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1994 lega->AddEntry(beforee,"Before correction","p");
1995 lega->AddEntry(aftere,"After correction","p");
1997 cV0Efficiencye->cd(1);
1999 TH1D *afteree = (TH1D *) aftere->Clone();
2000 afteree->SetMarkerStyle(stylee[binc]);
2001 afteree->SetMarkerColor(colorr[binc]);
2002 if(binc==0) afteree->Draw();
2003 else afteree->Draw("same");
2004 legtotal->AddEntry(afteree,(const char*) titlee,"p");
2005 cV0Efficiencye->cd(2);
2007 TH1D *aftereeu = (TH1D *) aftere->Clone();
2008 aftereeu->SetMarkerStyle(stylee[binc]);
2009 aftereeu->SetMarkerColor(colorr[binc]);
2010 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
2011 if(binc==0) aftereeu->Draw();
2012 else aftereeu->Draw("same");
2013 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
2014 ccV0Efficiency->cd(2);
2015 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
2016 efficiencyDDproj->SetTitle("");
2017 efficiencyDDproj->SetStats(0);
2018 efficiencyDDproj->SetMarkerStyle(25);
2019 efficiencyDDproj->Draw();
2022 cV0Efficiencye->cd(1);
2023 legtotal->Draw("same");
2024 cV0Efficiencye->cd(2);
2025 legtotalg->Draw("same");
2027 cenaxisa->SetRange(0,nbbin);
2028 cenaxisb->SetRange(0,nbbin);
2029 cenaxisc->SetRange(0,nbbin);
2031 if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
2040 //____________________________________________________________
2041 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
2044 // Unfold and eventually correct for efficiency the bgsubspectrum
2047 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2049 AliError("MC Container not available");
2054 AliError("No Correlation map available");
2059 AliCFDataGrid *dataGrid = 0x0;
2061 dataGrid = bgsubpectrum;
2065 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2067 AliError("Data Container not available");
2071 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2075 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2076 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2079 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2080 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2082 if(!fBeauty2ndMethod)
2084 // Consider parameterized IP cut efficiency
2085 if(!fInclusiveSpectrum){
2086 Int_t* bins=new Int_t[1];
2089 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2090 beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
2091 efficiencyD->Multiply(beffContainer,1);
2098 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2099 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2100 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2101 if(fSetSmoothing) unfolding.UseSmoothing();
2105 THnSparse* result = unfolding.GetUnfolded();
2106 THnSparse* residual = unfolding.GetEstMeasured();
2107 TList *listofresults = new TList;
2108 listofresults->SetOwner();
2109 listofresults->AddAt((THnSparse*)result->Clone(),0);
2110 listofresults->AddAt((THnSparse*)residual->Clone(),1);
2112 if(fDebugLevel > 0) {
2115 if(fBeamType==0) ptpr=0;
2116 if(fBeamType==1) ptpr=1;
2118 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
2119 cresidual->Divide(2,1);
2122 TGraphErrors* residualspectrumD = Normalize(residual);
2123 if(!residualspectrumD) {
2124 AliError("Number of Events not set for the normalization");
2127 residualspectrumD->SetTitle("");
2128 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2129 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2130 residualspectrumD->SetMarkerStyle(26);
2131 residualspectrumD->SetMarkerColor(kBlue);
2132 residualspectrumD->SetLineColor(kBlue);
2133 residualspectrumD->Draw("AP");
2134 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2135 dataGridBis->SetName("dataGridBis");
2136 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2137 measuredspectrumD->SetTitle("");
2138 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2139 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2140 measuredspectrumD->SetMarkerStyle(25);
2141 measuredspectrumD->SetMarkerColor(kBlack);
2142 measuredspectrumD->SetLineColor(kBlack);
2143 measuredspectrumD->Draw("P");
2144 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2145 legres->AddEntry(residualspectrumD,"Residual","p");
2146 legres->AddEntry(measuredspectrumD,"Measured","p");
2147 legres->Draw("same");
2149 TH1D *residualTH1D = residual->Projection(ptpr);
2150 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
2151 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2152 ratioresidual->SetName("ratioresidual");
2153 ratioresidual->SetTitle("");
2154 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2155 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2156 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2157 ratioresidual->SetStats(0);
2158 ratioresidual->Draw();
2160 if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
2163 return listofresults;
2167 //____________________________________________________________
2168 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
2171 // Apply unfolding and efficiency correction together to bgsubspectrum
2174 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
2176 AliError("MC Container not available");
2181 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2182 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2185 if(!fBeauty2ndMethod)
2187 // Consider parameterized IP cut efficiency
2188 if(!fInclusiveSpectrum){
2189 Int_t* bins=new Int_t[1];
2192 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2193 beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
2194 efficiencyD->Multiply(beffContainer,1);
2198 // Data in the right format
2199 AliCFDataGrid *dataGrid = 0x0;
2201 dataGrid = bgsubpectrum;
2205 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2207 AliError("Data Container not available");
2211 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2215 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2216 result->ApplyEffCorrection(*efficiencyD);
2218 if(fDebugLevel > 0) {
2221 if(fBeamType==0) ptpr=0;
2222 if(fBeamType==1) ptpr=1;
2224 printf("Step MC: %d\n",fStepMC);
2225 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2226 printf("Step MC true: %d\n",fStepTrue);
2227 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2228 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2229 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2231 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2233 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2235 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
2236 efficiencymcPIDD->SetTitle("");
2237 efficiencymcPIDD->SetStats(0);
2238 efficiencymcPIDD->SetMarkerStyle(25);
2239 efficiencymcPIDD->Draw();
2240 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
2241 efficiencymctrackinggeoD->SetTitle("");
2242 efficiencymctrackinggeoD->SetStats(0);
2243 efficiencymctrackinggeoD->SetMarkerStyle(26);
2244 efficiencymctrackinggeoD->Draw("same");
2245 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
2246 efficiencymcallD->SetTitle("");
2247 efficiencymcallD->SetStats(0);
2248 efficiencymcallD->SetMarkerStyle(27);
2249 efficiencymcallD->Draw("same");
2250 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
2251 efficiencyesdallD->SetTitle("");
2252 efficiencyesdallD->SetStats(0);
2253 efficiencyesdallD->SetMarkerStyle(24);
2254 efficiencyesdallD->Draw("same");
2255 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2256 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2257 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2258 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2259 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2260 legeff->Draw("same");
2264 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2265 TAxis *cenaxisa = sparseafter->GetAxis(0);
2266 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2267 TAxis *cenaxisb = sparsebefore->GetAxis(0);
2268 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2269 TAxis *cenaxisc = efficiencya->GetAxis(0);
2270 //Int_t nbbin = cenaxisb->GetNbins();
2271 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2272 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2273 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
2274 TString titlee("Efficiency_centrality_bin_");
2275 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2277 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2278 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2279 cefficiencye->Divide(2,1);
2280 cefficiencye->cd(1);
2282 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2283 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2284 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2285 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2286 CorrectFromTheWidth(afterefficiency);
2287 CorrectFromTheWidth(beforeefficiency);
2288 afterefficiency->SetStats(0);
2289 afterefficiency->SetTitle((const char*)titlee);
2290 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2291 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2292 afterefficiency->SetMarkerStyle(25);
2293 afterefficiency->SetMarkerColor(kBlack);
2294 afterefficiency->SetLineColor(kBlack);
2295 beforeefficiency->SetStats(0);
2296 beforeefficiency->SetTitle((const char*)titlee);
2297 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2298 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2299 beforeefficiency->SetMarkerStyle(24);
2300 beforeefficiency->SetMarkerColor(kBlue);
2301 beforeefficiency->SetLineColor(kBlue);
2302 afterefficiency->Draw();
2303 beforeefficiency->Draw("same");
2304 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2305 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2306 lega->AddEntry(afterefficiency,"After efficiency correction","p");
2308 cefficiencye->cd(2);
2309 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2310 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2311 efficiencyDDproj->SetTitle("");
2312 efficiencyDDproj->SetStats(0);
2313 efficiencyDDproj->SetMarkerStyle(25);
2314 efficiencyDDproj->SetMarkerColor(2);
2315 efficiencyDDproj->Draw();
2316 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2317 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2318 efficiencyDDproja->SetTitle("");
2319 efficiencyDDproja->SetStats(0);
2320 efficiencyDDproja->SetMarkerStyle(26);
2321 efficiencyDDproja->SetMarkerColor(4);
2322 efficiencyDDproja->Draw("same");
2326 if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
2333 //__________________________________________________________________________________
2334 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
2336 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2337 // Give the final pt spectrum to be compared
2340 if(fNEvents[i] > 0) {
2343 if(fBeamType==0) ptpr=0;
2344 if(fBeamType==1) ptpr=1;
2346 TH1D* projection = spectrum->Projection(ptpr);
2347 CorrectFromTheWidth(projection);
2348 TGraphErrors *graphError = NormalizeTH1(projection,i);
2357 //__________________________________________________________________________________
2358 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
2360 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2361 // Give the final pt spectrum to be compared
2363 if(fNEvents[i] > 0) {
2366 if(fBeamType==0) ptpr=0;
2367 if(fBeamType==1) ptpr=1;
2369 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
2370 CorrectFromTheWidth(projection);
2371 TGraphErrors *graphError = NormalizeTH1(projection,i);
2381 //__________________________________________________________________________________
2382 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2384 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2385 // Give the final pt spectrum to be compared
2387 Double_t chargecoefficient = 0.5;
2388 if(fChargeChoosen != 0) chargecoefficient = 1.0;
2390 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2391 printf("Normalizing Eta Range %f\n", etarange);
2392 if(fNEvents[i] > 0) {
2394 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2395 Double_t p = 0, dp = 0; Int_t point = 1;
2396 Double_t n = 0, dN = 0;
2397 Double_t nCorr = 0, dNcorr = 0;
2398 //Double_t errdN = 0, errdp = 0;
2400 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2401 point = ibin - input->GetXaxis()->GetFirst();
2402 p = input->GetXaxis()->GetBinCenter(ibin);
2403 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2404 n = input->GetBinContent(ibin);
2405 dN = input->GetBinError(ibin);
2407 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
2408 errdN = 1./(2. * TMath::Pi() * p);
2409 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2410 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2411 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN);
2413 spectrumNormalized->SetPoint(point, p, nCorr);
2414 spectrumNormalized->SetPointError(point, dp, dNcorr);
2416 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2417 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2418 spectrumNormalized->SetMarkerStyle(22);
2419 spectrumNormalized->SetMarkerColor(kBlue);
2420 spectrumNormalized->SetLineColor(kBlue);
2422 return spectrumNormalized;
2430 //__________________________________________________________________________________
2431 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
2433 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2434 // Give the final pt spectrum to be compared
2436 Double_t chargecoefficient = 0.5;
2437 if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
2439 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2440 printf("Normalizing Eta Range %f\n", etarange);
2441 if(normalization > 0) {
2443 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2444 Double_t p = 0, dp = 0; Int_t point = 1;
2445 Double_t n = 0, dN = 0;
2446 Double_t nCorr = 0, dNcorr = 0;
2447 //Double_t errdN = 0, errdp = 0;
2449 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2450 point = ibin - input->GetXaxis()->GetFirst();
2451 p = input->GetXaxis()->GetBinCenter(ibin);
2452 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2453 n = input->GetBinContent(ibin);
2454 dN = input->GetBinError(ibin);
2456 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
2457 errdN = 1./(2. * TMath::Pi() * p);
2458 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2459 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2460 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN);
2462 spectrumNormalized->SetPoint(point, p, nCorr);
2463 spectrumNormalized->SetPointError(point, dp, dNcorr);
2465 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2466 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2467 spectrumNormalized->SetMarkerStyle(22);
2468 spectrumNormalized->SetMarkerColor(kBlue);
2469 spectrumNormalized->SetLineColor(kBlue);
2471 return spectrumNormalized;
2479 //____________________________________________________________
2480 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2482 // Set the container for a given step to the
2484 if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
2485 fCFContainers->AddAt(cont, type);
2488 //____________________________________________________________
2489 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2491 // Get Correction Framework Container for given type
2493 if(!fCFContainers) return NULL;
2494 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2496 //____________________________________________________________
2497 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
2499 // Slice bin for a given source of electron
2500 // nDim is the number of dimension the corrections are done
2501 // dimensions are the definition of the dimensions
2502 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2503 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2504 // centrality (-1 means we do not cut on centrality)
2507 Double_t *varMin = new Double_t[container->GetNVar()],
2508 *varMax = new Double_t[container->GetNVar()];
2510 Double_t *binLimits;
2511 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2513 binLimits = new Double_t[container->GetNBins(ivar)+1];
2514 container->GetBinLimits(ivar,binLimits);
2515 varMin[ivar] = binLimits[0];
2516 varMax[ivar] = binLimits[container->GetNBins(ivar)];
2519 if((source>= 0) && (source<container->GetNBins(ivar))) {
2520 varMin[ivar] = container->GetAxis(4,0)->GetBinLowEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2521 varMax[ivar] = container->GetAxis(4,0)->GetBinUpEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2526 if(charge != kAllCharge){
2527 varMin[ivar] = container->GetAxis(3,0)->GetBinLowEdge(container->GetAxis(3,0)->FindBin(charge));
2528 varMax[ivar] = container->GetAxis(3,0)->GetBinUpEdge(container->GetAxis(3,0)->FindBin(charge));
2533 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2534 AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
2536 varMin[ivar] = fEtaRange[0];
2537 varMax[ivar] = fEtaRange[1];
2541 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2542 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2543 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2546 if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
2547 varMin[ivar] = binLimits[centralitylow];
2548 varMax[ivar] = binLimits[centralityhigh];
2550 TAxis *axistest = container->GetAxis(5,0);
2551 printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
2552 printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
2553 Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
2554 Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
2555 printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
2565 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2566 AddTemporaryObject(k);
2567 delete[] varMin; delete[] varMax;
2573 //_________________________________________________________________________
2574 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
2576 // Slice correlation
2579 Int_t ndimensions = correlationmatrix->GetNdimensions();
2580 //printf("Number of dimension %d correlation map\n",ndimensions);
2581 if(ndimensions < (2*nDim)) {
2582 AliError("Problem in the dimensions");
2586 // Cut in centrality is centrality > -1
2587 if((centralitylow >=0) && (centralityhigh >=0)) {
2589 TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
2590 TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
2592 Int_t bins0 = axiscentrality0->GetNbins();
2593 Int_t bins1 = axiscentrality1->GetNbins();
2595 printf("Number of centrality bins: %d and %d\n",bins0,bins1);
2596 if(bins0 != bins1) {
2597 AliError("Problem in the dimensions");
2601 if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
2602 axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
2603 axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
2605 Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
2606 Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
2607 Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
2608 Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
2609 printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
2610 printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
2612 TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
2613 ctestcorrelation->cd(1);
2614 TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
2615 projection->Draw("colz");
2622 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2623 //printf("Number of dimension %d container\n",ndimensionsContainer);
2625 Int_t *dim = new Int_t[nDim*2];
2626 for(Int_t iter=0; iter < nDim; iter++){
2627 dim[iter] = dimensions[iter];
2628 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2629 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2632 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2638 //___________________________________________________________________________
2639 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2641 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2644 TAxis *axis = h1->GetXaxis();
2645 Int_t nbinX = h1->GetNbinsX();
2647 for(Int_t i = 1; i <= nbinX; i++) {
2649 Double_t width = axis->GetBinWidth(i);
2650 Double_t content = h1->GetBinContent(i);
2651 Double_t error = h1->GetBinError(i);
2652 h1->SetBinContent(i,content/width);
2653 h1->SetBinError(i,error/width);
2658 //___________________________________________________________________________
2659 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2661 // Correct statistical error
2664 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2665 Int_t nbinX = h1->GetNbinsX();
2667 for(Long_t i = 1; i <= nbinX; i++) {
2669 Float_t content = h1->GetBinContent(i);
2671 Float_t error = TMath::Sqrt(content);
2672 backgroundGrid->SetElementError(bins, error);
2677 //____________________________________________________________
2678 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2680 // Emulate garbage collection on class level
2682 if(!fTemporaryObjects) {
2683 fTemporaryObjects= new TList;
2684 fTemporaryObjects->SetOwner();
2686 fTemporaryObjects->Add(o);
2689 //____________________________________________________________
2690 void AliHFEspectrum::ClearObject(TObject *o){
2692 // Do a safe deletion
2694 if(fTemporaryObjects){
2695 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2702 //_________________________________________________________________________
2703 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2704 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2707 //_________________________________________________________________________
2708 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2710 // Create efficiency grid and calculate efficiency
2713 TString name("eff");
2716 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2717 eff->CalculateEfficiency(step,step0);
2721 //____________________________________________________________________________
2722 THnSparse* AliHFEspectrum::GetCharmWeights(){
2725 // Measured D->e based weighting factors
2728 const Int_t nDimpp=1;
2729 Int_t nBinpp[nDimpp] = {35};
2730 Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
2731 const Int_t nDimPbPb=2;
2732 Int_t nBinPbPb[nDimPbPb] = {11,35};
2733 Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2735 Int_t looppt=nBinpp[0];
2738 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2739 fWeightCharm->SetBinEdges(0, ptbinning1);
2743 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2744 fWeightCharm->SetBinEdges(1, ptbinning1);
2745 fWeightCharm->SetBinEdges(0, kCentralityRange);
2746 loopcentr=nBinPbPb[0];
2749 // Weighting factor for pp
2750 Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
2752 // Weighting factor for PbPb (0-20%)
2753 //Double_t weight[35]={0.641897, 0.640472, 0.615228, 0.650469, 0.737762, 0.847867, 1.009317, 1.158594, 1.307482, 1.476973, 1.551131, 1.677131, 1.785478, 1.888933, 2.017957, 2.074757, 1.926700, 1.869495, 1.546558, 1.222873, 1.160313, 0.903375, 0.799642, 0.706244, 0.705449, 0.599947, 0.719570, 0.499422, 0.703978, 0.477452, 0.325057, 0.093391, 0.096675, 0.000000, 0.000000};
2755 // Weighting factor for PbPb (40-80%)
2756 //Double_t weight[35]={0.181953, 0.173248, 0.166799, 0.182558, 0.206581, 0.236955, 0.279390, 0.329129, 0.365260, 0.423059, 0.452057, 0.482726, 0.462627, 0.537770, 0.584663, 0.579452, 0.587194, 0.499498, 0.443299, 0.398596, 0.376695, 0.322331, 0.260890, 0.374834, 0.249114, 0.310330, 0.287326, 0.243174, 0.758945, 0.138867, 0.170576, 0.107797, 0.011390, 0.000000, 0.000000};
2760 Double_t contents[2];
2762 for(int icentr=0; icentr<loopcentr; icentr++)
2764 for(int i=0; i<looppt; i++){
2765 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2776 fWeightCharm->Fill(contents,weight[i]);
2780 Int_t nDimSparse = fWeightCharm->GetNdimensions();
2781 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2782 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2784 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2785 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2786 nBins *= binsvar[iVar];
2789 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2790 // loop that sets 0 error in each bin
2791 for (Long_t iBin=0; iBin<nBins; iBin++) {
2792 Long_t bintmp = iBin ;
2793 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2794 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2795 bintmp /= binsvar[iVar] ;
2797 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
2803 return fWeightCharm;
2806 //____________________________________________________________________________
2807 void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
2809 // TOF PID efficiencies
2811 if(fBeamType==0) ptpr=0;
2812 if(fBeamType==1) ptpr=1;
2815 const Int_t nCentralitybinning=11; //number of centrality bins
2818 loopcentr=nCentralitybinning;
2821 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2822 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2823 TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
2825 TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
2826 cefficiencyParamtof->cd();
2828 AliCFContainer *mccontainermcD = 0x0;
2829 AliCFContainer *mccontaineresdD = 0x0;
2830 TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2831 TH1D* efficiencyTOFPIDD[nCentralitybinning];
2832 TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2833 TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
2834 Int_t source = -1; //get parameterized TOF PID efficiencies
2836 for(int icentr=0; icentr<loopcentr; icentr++) {
2838 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2839 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2840 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2841 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2843 // mb sample for double check
2844 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2845 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2846 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2847 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2849 // mb sample with reconstructed variables
2850 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2851 else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2852 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2853 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2855 // mb sample with reconstructed variables
2856 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2857 else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2858 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2859 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2862 efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2863 efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
2864 efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2865 efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
2866 efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2867 efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
2868 efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2869 efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
2871 //fit (mc enhenced sample)
2872 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2873 efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2874 efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2875 fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
2877 //fit (esd enhenced sample)
2878 efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2879 efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2880 fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2884 // draw (for PbPb, only 1st bin)
2886 efficiencysigTOFPIDD[0]->SetTitle("");
2887 efficiencysigTOFPIDD[0]->SetStats(0);
2888 efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2889 efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2890 efficiencysigTOFPIDD[0]->SetLineColor(2);
2891 efficiencysigTOFPIDD[0]->Draw();
2894 efficiencyTOFPIDD[0]->SetTitle("");
2895 efficiencyTOFPIDD[0]->SetStats(0);
2896 efficiencyTOFPIDD[0]->SetMarkerStyle(24);
2897 efficiencyTOFPIDD[0]->SetMarkerColor(4);
2898 efficiencyTOFPIDD[0]->SetLineColor(4);
2899 efficiencyTOFPIDD[0]->Draw("same");
2902 efficiencysigesdTOFPIDD[0]->SetTitle("");
2903 efficiencysigesdTOFPIDD[0]->SetStats(0);
2904 efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2905 efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2906 efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2907 efficiencysigesdTOFPIDD[0]->Draw("same");
2910 efficiencyesdTOFPIDD[0]->SetTitle("");
2911 efficiencyesdTOFPIDD[0]->SetStats(0);
2912 efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2913 efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2914 efficiencyesdTOFPIDD[0]->SetLineColor(1);
2915 efficiencyesdTOFPIDD[0]->Draw("same");
2918 if(fEfficiencyTOFPIDD[0]){
2919 fEfficiencyTOFPIDD[0]->SetLineColor(2);
2920 fEfficiencyTOFPIDD[0]->Draw("same");
2923 if(fEfficiencyesdTOFPIDD[0]){
2924 fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2925 fEfficiencyesdTOFPIDD[0]->Draw("same");
2928 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2929 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2930 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2931 legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2932 legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2933 legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2934 legtofeff->Draw("same");
2937 TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
2938 cefficiencyParamIP->cd();
2939 gStyle->SetOptStat(0);
2941 // IP cut efficiencies
2942 for(int icentr=0; icentr<loopcentr; icentr++) {
2944 AliCFContainer *charmCombined = 0x0;
2945 AliCFContainer *beautyCombined = 0x0;
2946 AliCFContainer *beautyCombinedesd = 0x0;
2948 printf("centrality printing %i \n",icentr);
2950 source = 0; //charm enhenced
2951 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2952 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2953 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2954 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2956 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
2958 source = 1; //beauty enhenced
2959 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2960 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2961 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2962 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2964 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
2966 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2967 else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2968 AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
2969 efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2971 beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
2973 source = 0; //charm mb
2974 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2975 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2976 AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2977 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2979 charmCombined->Add(mccontainermcD);
2980 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
2981 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2983 source = 1; //beauty mb
2984 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2985 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2986 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
2987 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2989 beautyCombined->Add(mccontainermcD);
2990 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
2991 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2993 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2994 else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2995 AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
2996 efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2998 beautyCombinedesd->Add(mccontaineresdD);
2999 AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
3000 efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3002 source = 2; //conversion mb
3003 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3004 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3005 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
3006 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3008 source = 3; //non HFE except for the conversion mb
3009 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3010 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3011 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
3012 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3014 if(fIPEffCombinedSamples){
3015 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
3016 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
3017 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
3020 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
3021 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
3022 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
3024 fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
3025 fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
3026 fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
3027 fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
3032 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
3033 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
3035 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
3036 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
3039 for(int icentr=0; icentr<loopcentr; icentr++) {
3040 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3041 fipfit->SetLineColor(2);
3042 fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
3043 fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3044 if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
3045 else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
3047 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3048 fipfit->SetLineColor(6);
3049 fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
3050 fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
3051 if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
3052 else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
3054 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3055 fipfit->SetLineColor(1);
3056 fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
3057 fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3058 fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
3060 if(fIPParameterizedEff){
3061 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3062 fipfitnonhfe->SetLineColor(3);
3063 fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
3064 fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3065 fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
3067 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3068 fipfitnonhfe->SetLineColor(4);
3069 fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
3070 fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3071 fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
3075 // draw (for PbPb, only 1st bin)
3076 fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
3077 fEfficiencyCharmSigD[0]->SetMarkerColor(1);
3078 fEfficiencyCharmSigD[0]->SetLineColor(1);
3079 fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
3080 fEfficiencyBeautySigD[0]->SetMarkerColor(2);
3081 fEfficiencyBeautySigD[0]->SetLineColor(2);
3082 fEfficiencyBeautySigesdD[0]->SetStats(0);
3083 fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
3084 fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
3085 fEfficiencyBeautySigesdD[0]->SetLineColor(6);
3086 fCharmEff[0]->SetMarkerStyle(24);
3087 fCharmEff[0]->SetMarkerColor(1);
3088 fCharmEff[0]->SetLineColor(1);
3089 fBeautyEff[0]->SetMarkerStyle(24);
3090 fBeautyEff[0]->SetMarkerColor(2);
3091 fBeautyEff[0]->SetLineColor(2);
3092 fConversionEff[0]->SetMarkerStyle(24);
3093 fConversionEff[0]->SetMarkerColor(3);
3094 fConversionEff[0]->SetLineColor(3);
3095 fNonHFEEff[0]->SetMarkerStyle(24);
3096 fNonHFEEff[0]->SetMarkerColor(4);
3097 fNonHFEEff[0]->SetLineColor(4);
3099 fEfficiencyCharmSigD[0]->Draw();
3100 fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
3101 fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
3103 fEfficiencyBeautySigD[0]->Draw("same");
3104 fEfficiencyBeautySigesdD[0]->Draw("same");
3105 //fCharmEff[0]->Draw("same");
3106 //fBeautyEff[0]->Draw("same");
3109 fConversionEffbgc->SetMarkerStyle(25);
3110 fConversionEffbgc->SetMarkerColor(3);
3111 fConversionEffbgc->SetLineColor(3);
3112 fNonHFEEffbgc->SetMarkerStyle(25);
3113 fNonHFEEffbgc->SetMarkerColor(4);
3114 fNonHFEEffbgc->SetLineColor(4);
3115 fConversionEffbgc->Draw("same");
3116 fNonHFEEffbgc->Draw("same");
3119 fConversionEff[0]->Draw("same");
3120 fNonHFEEff[0]->Draw("same");
3122 if(fEfficiencyIPBeautyD[0])
3123 fEfficiencyIPBeautyD[0]->Draw("same");
3124 if(fEfficiencyIPBeautyesdD[0])
3125 fEfficiencyIPBeautyesdD[0]->Draw("same");
3126 if( fEfficiencyIPCharmD[0])
3127 fEfficiencyIPCharmD[0]->Draw("same");
3128 if(fIPParameterizedEff){
3129 fEfficiencyIPConversionD[0]->Draw("same");
3130 fEfficiencyIPNonhfeD[0]->Draw("same");
3132 TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
3133 legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
3134 legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
3135 legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
3136 legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
3137 legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
3138 legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
3139 //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
3140 //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
3141 legipeff->Draw("same");
3143 //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
3146 //____________________________________________________________________________
3147 THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
3149 // Return beauty electron IP cut efficiency
3152 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3153 const Int_t nCentralitybinning=11;//number of centrality bins
3154 Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3155 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3157 Int_t nDim=1; //dimensions of the efficiency weighting grid
3161 nDim=2; //dimensions of the efficiency weighting grid
3163 Int_t nBin[1] = {nPtbinning1};
3164 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3168 if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
3169 else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3171 for(Int_t idim = 0; idim < nDim; idim++)
3173 if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
3176 ipcut->SetBinEdges(0, kCentralityRange);
3177 ipcut->SetBinEdges(1, kPtRange);
3181 Double_t weight[nCentralitybinning];
3182 Double_t weightErr[nCentralitybinning];
3183 Double_t contents[2];
3185 for(Int_t a=0;a<11;a++)
3192 Int_t looppt=nBin[0];
3197 loopcentr=nBinPbPb[0];
3201 for(int icentr=0; icentr<loopcentr; icentr++)
3203 for(int i=0; i<looppt; i++)
3205 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3207 if(fEfficiencyIPBeautyD[icentr]){
3208 weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
3209 weightErr[icentr] = 0;
3212 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3213 weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
3214 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3218 if(fEfficiencyIPBeautyesdD[icentr]){
3219 weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
3220 weightErr[icentr] = 0;
3223 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3224 weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
3225 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3239 ipcut->Fill(contents,weight[icentr]);
3240 ipcut->SetBinError(ibin,weightErr[icentr]);
3244 Int_t nDimSparse = ipcut->GetNdimensions();
3245 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3246 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3248 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3249 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
3250 nBins *= binsvar[iVar];
3253 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3254 // loop that sets 0 error in each bin
3255 for (Long_t iBin=0; iBin<nBins; iBin++) {
3256 Long_t bintmp = iBin ;
3257 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3258 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3259 bintmp /= binsvar[iVar] ;
3261 //ipcut->SetBinError(binfill,0.); // put 0 everywhere
3270 //____________________________________________________________________________
3271 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3273 // Return PID x IP cut efficiency
3275 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3276 const Int_t nCentralitybinning=11;//number of centrality bins
3277 Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3278 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3280 Int_t nDim=1; //dimensions of the efficiency weighting grid
3284 nDim=2; //dimensions of the efficiency weighting grid
3286 Int_t nBin[1] = {nPtbinning1};
3287 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3290 if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3291 else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3292 for(Int_t idim = 0; idim < nDim; idim++)
3295 if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3298 pideff->SetBinEdges(0, kCentralityRange);
3299 pideff->SetBinEdges(1, kPtRange);
3304 Double_t weight[nCentralitybinning];
3305 Double_t weightErr[nCentralitybinning];
3306 Double_t contents[2];
3308 for(Int_t a=0;a<nCentralitybinning;a++)
3314 Int_t looppt=nBin[0];
3319 loopcentr=nBinPbPb[0];
3322 for(int icentr=0; icentr<loopcentr; icentr++)
3324 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3325 for(int i=0; i<looppt; i++)
3327 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3329 Double_t tofpideff = 0.;
3330 Double_t tofpideffesd = 0.;
3331 if(fEfficiencyTOFPIDD[icentr])
3332 tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
3334 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3336 if(fEfficiencyesdTOFPIDD[icentr])
3337 tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3339 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3342 //tof pid eff x tpc pid eff x ip cut eff
3343 if(fIPParameterizedEff){
3345 if(fEfficiencyIPCharmD[icentr]){
3346 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3347 weightErr[icentr] = 0;
3350 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3351 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3352 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3355 else if(source==2) {
3356 if(fEfficiencyIPConversionD[icentr]){
3357 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
3358 weightErr[icentr] = 0;
3361 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
3362 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3363 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3366 else if(source==3) {
3367 if(fEfficiencyIPNonhfeD[icentr]){
3368 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
3369 weightErr[icentr] = 0;
3372 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
3373 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3374 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3380 if(fEfficiencyIPCharmD[icentr]){
3381 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3382 weightErr[icentr] = 0;
3385 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3386 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3387 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3392 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3393 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3396 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3397 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3402 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3403 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3406 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3407 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3423 pideff->Fill(contents,weight[icentr]);
3424 pideff->SetBinError(ibin,weightErr[icentr]);
3428 Int_t nDimSparse = pideff->GetNdimensions();
3429 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3430 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3432 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3433 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3434 nBins *= binsvar[iVar];
3437 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3438 // loop that sets 0 error in each bin
3439 for (Long_t iBin=0; iBin<nBins; iBin++) {
3440 Long_t bintmp = iBin ;
3441 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3442 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3443 bintmp /= binsvar[iVar] ;
3454 //__________________________________________________________________________
3455 AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
3457 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
3471 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
3472 const Int_t nCentralitybinning=11;//number of centrality bins
3473 Double_t kPtRange[19] = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 12, 16, 20};
3475 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3476 Int_t nBin[1] = {nPtbinning1};
3477 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3479 AliCFDataGrid *rawBeautyContainer;
3480 if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
3481 else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
3482 // printf("number of bins= %d\n",bins[0]);
3487 THnSparseF *brawspectra;
3488 if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
3489 else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
3490 if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
3493 // brawspectra->SetBinEdges(0, centralityBins);
3494 brawspectra->SetBinEdges(0, kCentralityRange);
3495 brawspectra->SetBinEdges(1, kPtRange);
3499 Double_t yields= 0.;
3500 Double_t valuesb[2];
3502 //Int_t looppt=nBin[0];
3506 loopcentr=nBinPbPb[0];
3509 for(int icentr=0; icentr<loopcentr; icentr++)
3512 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
3513 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3515 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
3522 if(fBeamType==0) valuesb[0]=pt[0];
3523 brawspectra->Fill(valuesb,yields);
3529 Int_t nDimSparse = brawspectra->GetNdimensions();
3530 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3531 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3533 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3534 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
3535 nBins *= binsvar[iVar];
3538 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3539 // loop that sets 0 error in each bin
3540 for (Long_t iBin=0; iBin<nBins; iBin++) {
3541 Long_t bintmp = iBin ;
3542 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3543 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3544 bintmp /= binsvar[iVar] ;
3546 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
3550 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
3551 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
3554 fBSpectrum2ndMethod->SetMarkerStyle(24);
3555 fBSpectrum2ndMethod->Draw("p");
3556 hRawBeautySpectra->SetMarkerStyle(25);
3557 hRawBeautySpectra->Draw("samep");
3562 return rawBeautyContainer;
3565 //__________________________________________________________________________
3566 void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
3568 // Calculate non HFE sys
3575 Double_t evtnorm[1] = {0.0};
3576 if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
3578 AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
3579 AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
3581 AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
3582 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
3583 AliCFDataGrid *bgConvGrid[kBgLevels];
3585 Int_t stepbackground = 3;
3586 Int_t* bins=new Int_t[1];
3587 const Char_t *bgBase[2] = {"pi0","eta"};
3589 bins[0]=fConversionEff[centrality]->GetNbinsX();
3591 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
3592 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
3594 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
3595 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
3596 convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground);
3597 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
3598 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
3600 nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground);
3601 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
3602 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
3605 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
3606 for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
3607 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
3610 bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
3612 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
3613 for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
3614 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
3617 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
3619 bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
3620 bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
3622 bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
3623 bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
3628 //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0)
3629 AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
3630 TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
3631 for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
3632 bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
3633 bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
3634 bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
3635 bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
3637 //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
3639 hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
3640 hBaseErrors[iErr][0]->Scale(-1.);
3641 hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
3642 hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
3647 //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
3648 TH1D *hSpeciesErrors[kElecBgSources-1];
3649 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3650 if(fEtaSyst && (iSource == 1))continue;
3651 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
3652 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
3653 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
3654 hSpeciesErrors[iSource-1]->Scale(0.3);
3657 //Int_t firstBgSource = 0;//if eta systematics are not from scaling
3658 //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
3659 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
3660 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
3661 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
3663 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
3664 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
3666 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
3667 Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
3668 pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
3669 pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
3671 etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
3672 etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
3674 else{ etaErrLow = etaErrUp = 0.;}
3676 Double_t sqrsumErrs= 0;
3677 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3678 if(fEtaSyst && (iSource == 1))continue;
3679 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
3680 sqrsumErrs+=(scalingErr*scalingErr);
3682 for(Int_t iErr = 0; iErr < 2; iErr++){
3683 for(Int_t iLevel = 0; iLevel < 2; iLevel++){
3684 hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
3688 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
3689 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
3690 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
3692 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
3693 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
3697 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
3699 cPiErrors->SetLogx();
3700 cPiErrors->SetLogy();
3701 hBaseErrors[0][0]->Draw();
3702 //hBaseErrors[0][1]->SetMarkerColor(kBlack);
3703 //hBaseErrors[0][1]->SetLineColor(kBlack);
3704 //hBaseErrors[0][1]->Draw("SAME");
3706 hBaseErrors[1][0]->Draw("SAME");
3707 hBaseErrors[1][0]->SetMarkerColor(kBlack);
3708 hBaseErrors[1][0]->SetLineColor(kBlack);
3709 //hBaseErrors[1][1]->SetMarkerColor(13);
3710 //hBaseErrors[1][1]->SetLineColor(13);
3711 //hBaseErrors[1][1]->Draw("SAME");
3713 //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
3714 //hOverallBinScaledErrsUp->SetLineColor(kBlue);
3715 //hOverallBinScaledErrsUp->Draw("SAME");
3716 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
3717 hOverallBinScaledErrsLow->SetLineColor(kGreen);
3718 hOverallBinScaledErrsLow->Draw("SAME");
3719 hScalingErrors->SetLineColor(kBlue);
3720 hScalingErrors->Draw("SAME");
3722 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
3723 lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
3724 //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
3726 lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
3727 //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
3729 lPiErr->AddEntry(hScalingErrors, "scaling error");
3730 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
3731 //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
3732 lPiErr->Draw("SAME");
3735 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
3736 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
3737 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
3738 hLowSystScaled->Scale(evtnorm[0]);
3739 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
3740 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
3741 //histograms to be normalized to TGraphErrors
3742 CorrectFromTheWidth(hNormAllSystErrUp);
3743 CorrectFromTheWidth(hNormAllSystErrLow);
3745 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
3747 cNormOvErrs->SetLogx();
3748 cNormOvErrs->SetLogy();
3750 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
3751 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
3752 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
3753 gOverallSystErrUp->SetMarkerColor(kBlack);
3754 gOverallSystErrUp->SetLineColor(kBlack);
3755 gOverallSystErrLow->SetMarkerColor(kRed);
3756 gOverallSystErrLow->SetLineColor(kRed);
3757 gOverallSystErrUp->Draw("AP");
3758 gOverallSystErrLow->Draw("PSAME");
3759 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
3760 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
3761 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
3762 lAllSys->Draw("same");
3765 AliCFDataGrid *bgYieldGrid;
3767 bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
3769 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
3771 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
3772 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
3773 hRelErrUp->Divide(hBgYield);
3774 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
3775 hRelErrLow->Divide(hBgYield);
3777 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
3779 cRelErrs->SetLogx();
3780 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
3782 hRelErrLow->SetLineColor(kBlack);
3783 hRelErrLow->Draw("SAME");
3785 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
3786 lRel->AddEntry(hRelErrUp, "upper");
3787 lRel->AddEntry(hRelErrLow, "lower");
3790 //CorrectFromTheWidth(hBgYield);
3791 //hBgYield->Scale(evtnorm[0]);
3794 //write histograms/TGraphs to file
3795 TFile *output = new TFile("systHists.root","recreate");
3797 hBgYield->SetName("hBgYield");
3799 hRelErrUp->SetName("hRelErrUp");
3801 hRelErrLow->SetName("hRelErrLow");
3802 hRelErrLow->Write();
3803 hUpSystScaled->SetName("hOverallSystErrUp");
3804 hUpSystScaled->Write();
3805 hLowSystScaled->SetName("hOverallSystErrLow");
3806 hLowSystScaled->Write();
3807 gOverallSystErrUp->SetName("gOverallSystErrUp");
3808 gOverallSystErrUp->Write();
3809 gOverallSystErrLow->SetName("gOverallSystErrLow");
3810 gOverallSystErrLow->Write();
3817 //____________________________________________________________
3818 void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
3821 // Unfold backgrounds to check its sanity
3824 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
3825 //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
3827 AliError("MC Container not available");
3832 AliError("No Correlation map available");
3837 AliCFDataGrid *dataGrid = 0x0;
3838 dataGrid = bgsubpectrum;
3841 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
3842 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
3845 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
3846 efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
3848 // Unfold background spectra
3850 if(fBeamType==0)nDim = 1;
3851 if(fBeamType==1)nDim = 2;
3852 AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
3853 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
3854 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
3855 if(fSetSmoothing) unfolding.UseSmoothing();
3859 THnSparse* result = unfolding.GetUnfolded();
3860 TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
3865 result->GetAxis(0)->SetRange(1,1);
3866 TH1D* htest1=(TH1D*)result->Projection(0);
3869 result->GetAxis(0)->SetRange(1,1);
3870 TH1D* htest2=(TH1D*)result->Projection(1);
3873 result->GetAxis(0)->SetRange(6,6);
3874 TH1D* htest3=(TH1D*)result->Projection(0);
3877 result->GetAxis(0)->SetRange(6,6);
3878 TH1D* htest4=(TH1D*)result->Projection(1);
3887 TGraphErrors* unfoldedbgspectrumD = Normalize(result);
3888 if(!unfoldedbgspectrumD) {
3889 AliError("Unfolded background spectrum doesn't exist");
3892 TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
3893 if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
3898 result->GetAxis(0)->SetRange(centr,centr);
3899 unfoldedbgspectrumD = Normalize(result,centr-1);
3900 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
3902 result->GetAxis(0)->SetRange(centr,centr);
3903 unfoldedbgspectrumD = Normalize(result,centr-1);
3904 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");