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),
94 fChargeChoosen(kAllCharge),
95 fNCentralityBinAtTheEnd(0),
96 fTestCentralityLow(-1),
97 fTestCentralityHigh(-1),
98 fFillMoreCorrelationMatrix(kFALSE),
99 fHadronEffbyIPcut(NULL),
100 fConversionEffbgc(NULL),
102 fBSpectrum2ndMethod(NULL),
103 fkBeauty2ndMethodfilename(""),
107 fWriteToFile(kFALSE),
111 // Default constructor
114 for(Int_t k = 0; k < 20; k++){
117 fLowBoundaryCentralityBinAtTheEnd[k] = 0;
118 fHighBoundaryCentralityBinAtTheEnd[k] = 0;
121 fEfficiencyTOFPIDD[k] = 0;
122 fEfficiencyesdTOFPIDD[k] = 0;
123 fEfficiencyIPCharmD[k] = 0;
124 fEfficiencyIPBeautyD[k] = 0;
125 fEfficiencyIPBeautyesdD[k] = 0;
126 fEfficiencyIPConversionD[k] = 0;
127 fEfficiencyIPNonhfeD[k] = 0;
129 fConversionEff[k] = 0;
133 fEfficiencyCharmSigD[k] = 0;
134 fEfficiencyBeautySigD[k] = 0;
135 fEfficiencyBeautySigesdD[k] = 0;
138 memset(fEtaRange, 0, sizeof(Double_t) * 2);
139 memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
140 memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
141 memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
144 //____________________________________________________________
145 AliHFEspectrum::~AliHFEspectrum(){
149 if(fCFContainers) delete fCFContainers;
150 if(fTemporaryObjects){
151 fTemporaryObjects->Clear();
152 delete fTemporaryObjects;
155 //____________________________________________________________
156 Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
158 // Init what we need for the correction:
160 // Raw spectrum, hadron contamination
161 // MC efficiency maps, correlation matrix
162 // V0 efficiency if wanted
164 // This for a given dimension.
165 // If no inclusive spectrum, then take only efficiency map for beauty electron
166 // and the appropriate correlation matrix
170 if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
175 if(fBeamType==0) kNdim=3;
184 // Get the requested format
187 switch(fNbDimensions){
190 case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
192 case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
195 AliError("Container with this number of dimensions not foreseen (yet)");
201 // PbPb analysis; centrality as first dimension
202 Int_t nbDimensions = fNbDimensions;
203 fNbDimensions = fNbDimensions + 1;
204 switch(nbDimensions){
209 for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
212 for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
215 AliError("Container with this number of dimensions not foreseen (yet)");
220 // Data container: raw spectrum + hadron contamination
221 AliCFContainer *datacontainer = 0x0;
222 if(fInclusiveSpectrum) {
223 datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
227 datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
229 AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
230 if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
232 AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
233 AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
234 if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
236 SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
237 SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
239 // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
240 AliCFContainer *mccontaineresd = 0x0;
241 AliCFContainer *mccontaineresdbg = 0x0;
242 AliCFContainer *mccontainermc = 0x0;
243 AliCFContainer *mccontainermcbg = 0x0;
244 AliCFContainer *nonHFEweightedContainer = 0x0;
245 AliCFContainer *convweightedContainer = 0x0;
246 AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
247 AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
249 if(fInclusiveSpectrum) {
250 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
251 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
254 mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
255 mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
256 mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
257 mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
260 const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
261 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
262 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
263 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
264 nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
265 convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
266 for(Int_t icentr=0;icentr<kNcentr;icentr++)
270 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
271 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
276 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
277 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
279 // if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
281 if(fBeamType == 1)break;
286 nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
287 convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
288 if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
291 if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
294 if(!fInclusiveSpectrum) source = 1; //beauty
295 AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
296 AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
297 if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
298 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
299 SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
301 // set charm, nonHFE container to estimate BG
302 if(!fInclusiveSpectrum){
304 mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
305 SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
308 AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
309 SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
310 AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
311 SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
314 SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
317 // MC container: correlation matrix
318 THnSparseF *mccorrelation = 0x0;
319 if(fInclusiveSpectrum) {
320 if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
321 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
322 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
323 else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
324 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
326 if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
328 else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
329 //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
330 if(!mccorrelation) return kFALSE;
331 Int_t testCentralityLow = fTestCentralityLow;
332 Int_t testCentralityHigh = fTestCentralityHigh;
333 if(fFillMoreCorrelationMatrix) {
334 testCentralityLow = fTestCentralityLow-1;
335 testCentralityHigh = fTestCentralityHigh+1;
337 THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
338 if(!mccorrelationD) {
339 printf("No correlation\n");
342 SetCorrelation(mccorrelationD);
344 // V0 container Electron, pt eta phi
346 AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
347 if(!containerV0) return kFALSE;
348 AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
349 if(!containerV0Electron) return kFALSE;
350 SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
356 AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
357 contaminationspectrum->SetName("contaminationspectrum");
358 TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
359 ccontaminationspectrum->Divide(3,1);
360 ccontaminationspectrum->cd(1);
361 TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
362 TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
363 TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
364 contaminationspectrum2dpteta->SetStats(0);
365 contaminationspectrum2dpteta->SetTitle("");
366 contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
367 contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
368 contaminationspectrum2dptphi->SetStats(0);
369 contaminationspectrum2dptphi->SetTitle("");
370 contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
371 contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
372 contaminationspectrum2detaphi->SetStats(0);
373 contaminationspectrum2detaphi->SetTitle("");
374 contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
375 contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
376 contaminationspectrum2dptphi->Draw("colz");
377 ccontaminationspectrum->cd(2);
378 contaminationspectrum2dpteta->Draw("colz");
379 ccontaminationspectrum->cd(3);
380 contaminationspectrum2detaphi->Draw("colz");
382 TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
386 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
387 ptcorrelation->Draw("colz");
390 TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
391 ptcorrelation->Draw("colz");
395 ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
396 ccorrelation->SaveAs("correlationMatrix.eps");
400 TFile *file = TFile::Open("tests.root","recreate");
401 datacontainerD->Write("data");
402 mccontainermcD->Write("mcefficiency");
403 mccorrelationD->Write("correlationmatrix");
410 //____________________________________________________________
411 void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
413 // get spectrum for beauty 2nd method
416 TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
417 fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
421 //____________________________________________________________
422 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
424 // Correct the spectrum for efficiency and unfolding
425 // with both method and compare
428 gStyle->SetPalette(1);
429 gStyle->SetOptStat(1111);
430 gStyle->SetPadBorderMode(0);
431 gStyle->SetCanvasColor(10);
432 gStyle->SetPadLeftMargin(0.13);
433 gStyle->SetPadRightMargin(0.13);
435 printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
437 ///////////////////////////
438 // Check initialization
439 ///////////////////////////
441 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
442 AliInfo("You have to init before");
446 if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
447 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
451 SetNumberOfIteration(10);
452 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
454 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
455 //////////////////////////////////
456 // Subtract hadron background
457 /////////////////////////////////
458 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
459 if(subtractcontamination) {
460 dataspectrumaftersubstraction = SubtractBackground(kTRUE);
461 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
464 printf("cloning spectrum\n");
465 AliCFDataGrid *rawsave(NULL);
466 if(dataspectrumaftersubstraction){
467 rawsave = (AliCFDataGrid *)dataspectrumaftersubstraction->Clone("rawdata");
469 AliCFContainer *dataContainer = GetContainer(kDataContainer);
471 AliError("Data Container not available");
473 rawsave = new AliCFDataGrid("rawsave", "raw spectrum after subtraction",*dataContainer, fStepData);
475 printf("cloned: %p\n", rawsave);
477 ////////////////////////////////////////////////
478 // Correct for TPC efficiency from V0
479 ///////////////////////////////////////////////
480 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
481 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
482 if(dataContainerV0){printf("Got the V0 container\n");
483 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
484 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
487 //////////////////////////////////////////////////////////////////////////////
488 // Correct for efficiency parametrized (if TPC efficiency is parametrized)
489 /////////////////////////////////////////////////////////////////////////////
490 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
491 if(fEfficiencyFunction){
492 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
493 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
499 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
501 printf("Unfolded failed\n");
504 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
505 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
506 if(!correctedspectrum){
507 AliError("No corrected spectrum\n");
510 if(!residualspectrum){
511 AliError("No residul spectrum\n");
515 /////////////////////
518 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
524 if(fDebugLevel > 0.0) {
527 if(fBeamType==0) ptpr=0;
528 if(fBeamType==1) ptpr=1;
530 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
531 ccorrected->Divide(2,1);
534 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
535 correctedspectrumD->SetTitle("");
536 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
537 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
538 correctedspectrumD->SetMarkerStyle(26);
539 correctedspectrumD->SetMarkerColor(kBlue);
540 correctedspectrumD->SetLineColor(kBlue);
541 correctedspectrumD->Draw("AP");
542 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
543 alltogetherspectrumD->SetTitle("");
544 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
545 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
546 alltogetherspectrumD->SetMarkerStyle(25);
547 alltogetherspectrumD->SetMarkerColor(kBlack);
548 alltogetherspectrumD->SetLineColor(kBlack);
549 alltogetherspectrumD->Draw("P");
550 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
551 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
552 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
553 legcorrected->Draw("same");
555 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
556 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
557 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
558 ratiocorrected->SetName("ratiocorrected");
559 ratiocorrected->SetTitle("");
560 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
561 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
562 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
563 ratiocorrected->SetStats(0);
564 ratiocorrected->Draw();
565 if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
567 //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
568 //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
569 //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
570 //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
572 TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
573 TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
574 TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
575 TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
579 TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
580 ccorrectedallspectra->Divide(2,1);
581 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
582 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
584 THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
585 TAxis *cenaxisa = sparsesu->GetAxis(0);
586 THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
587 TAxis *cenaxisb = sparsed->GetAxis(0);
588 Int_t nbbin = cenaxisb->GetNbins();
589 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
590 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
591 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
592 TString titlee("corrected_centrality_bin_");
594 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
596 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
598 TString titleec("corrected_check_projection_bin_");
600 titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
602 titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
604 TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
605 titleenameunotnormalized += "[";
606 titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
607 titleenameunotnormalized += "_";
608 titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
609 titleenameunotnormalized += "[";
610 TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
611 titleenameunormalized += "[";
612 titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
613 titleenameunormalized += "_";
614 titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
615 titleenameunormalized += "[";
616 TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
617 titleenamednotnormalized += "[";
618 titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
619 titleenamednotnormalized += "_";
620 titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
621 titleenamednotnormalized += "[";
622 TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
623 titleenamednormalized += "[";
624 titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
625 titleenamednormalized += "_";
626 titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
627 titleenamednormalized += "[";
629 for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
630 printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
631 nbEvents += fNEvents[k];
633 Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
634 Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
635 printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
636 Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
637 Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
638 printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
639 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
640 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
641 TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
642 ccorrectedcheck->cd(1);
643 TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
644 TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
646 aftersdc->Draw("same");
647 TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
648 ccorrectede->Divide(2,1);
651 TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
652 CorrectFromTheWidth(aftersu);
653 aftersu->SetName((const char*)titleenameunotnormalized);
654 unfoldingspectrac[binc] = *aftersu;
656 TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
657 aftersun->SetTitle("");
658 aftersun->GetYaxis()->SetTitleOffset(1.5);
659 aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
660 aftersun->SetMarkerStyle(26);
661 aftersun->SetMarkerColor(kBlue);
662 aftersun->SetLineColor(kBlue);
663 aftersun->Draw("AP");
664 aftersun->SetName((const char*)titleenameunormalized);
665 unfoldingspectracn[binc] = *aftersun;
667 TH1D *aftersd = (TH1D *) sparsed->Projection(1);
668 CorrectFromTheWidth(aftersd);
669 aftersd->SetName((const char*)titleenamednotnormalized);
670 correctedspectrac[binc] = *aftersd;
672 TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
673 aftersdn->SetTitle("");
674 aftersdn->GetYaxis()->SetTitleOffset(1.5);
675 aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
676 aftersdn->SetMarkerStyle(25);
677 aftersdn->SetMarkerColor(kBlack);
678 aftersdn->SetLineColor(kBlack);
680 aftersdn->SetName((const char*)titleenamednormalized);
681 correctedspectracn[binc] = *aftersdn;
682 TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
683 legcorrectedud->AddEntry(aftersun,"Corrected","p");
684 legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
685 legcorrectedud->Draw("same");
686 ccorrectedallspectra->cd(1);
688 TH1D *aftersunn = (TH1D *) aftersun->Clone();
689 aftersunn->SetMarkerStyle(stylee[binc]);
690 aftersunn->SetMarkerColor(colorr[binc]);
691 if(binc==0) aftersunn->Draw("AP");
692 else aftersunn->Draw("P");
693 legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
694 ccorrectedallspectra->cd(2);
696 TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
697 aftersdnn->SetMarkerStyle(stylee[binc]);
698 aftersdnn->SetMarkerColor(colorr[binc]);
699 if(binc==0) aftersdnn->Draw("AP");
700 else aftersdnn->Draw("P");
701 legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
703 TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
704 TString titleee("ratiocorrected_bin_");
706 ratiocorrectedbinc->SetName((const char*) titleee);
707 ratiocorrectedbinc->SetTitle("");
708 ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
709 ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
710 ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
711 ratiocorrectedbinc->SetStats(0);
712 ratiocorrectedbinc->Draw();
715 ccorrectedallspectra->cd(1);
716 legtotal->Draw("same");
717 ccorrectedallspectra->cd(2);
718 legtotalg->Draw("same");
720 cenaxisa->SetRange(0,nbbin);
721 cenaxisb->SetRange(0,nbbin);
722 if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
725 // Dump to file if needed
727 TFile *out = new TFile("finalSpectrum.root","recreate");
728 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
729 correctedspectrumD->Write();
730 alltogetherspectrumD->SetName("AlltogetherSpectrum");
731 alltogetherspectrumD->Write();
732 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
733 ratiocorrected->Write();
734 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
735 correctedspectrum->Write();
736 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
737 alltogetherCorrection->Write();
739 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
740 unfoldingspectrac[binc].Write();
741 unfoldingspectracn[binc].Write();
742 correctedspectrac[binc].Write();
743 correctedspectracn[binc].Write();
745 out->Close(); delete out;
748 if (unfoldingspectrac) delete[] unfoldingspectrac;
749 if (unfoldingspectracn) delete[] unfoldingspectracn;
750 if (correctedspectrac) delete[] correctedspectrac;
751 if (correctedspectracn) delete[] correctedspectracn;
758 //____________________________________________________________
759 Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
761 // Correct the spectrum for efficiency and unfolding for beauty analysis
762 // with both method and compare
765 gStyle->SetPalette(1);
766 gStyle->SetOptStat(1111);
767 gStyle->SetPadBorderMode(0);
768 gStyle->SetCanvasColor(10);
769 gStyle->SetPadLeftMargin(0.13);
770 gStyle->SetPadRightMargin(0.13);
772 ///////////////////////////
773 // Check initialization
774 ///////////////////////////
776 if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
777 AliInfo("You have to init before");
781 if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
782 AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
786 SetNumberOfIteration(10);
787 SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
789 AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
790 //////////////////////////////////
791 // Subtract hadron background
792 /////////////////////////////////
793 AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
794 AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
795 TGraphErrors *gNormalizedRawSpectrum = 0x0;
796 if(subtractcontamination) {
797 if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
798 else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
799 unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
800 dataGridAfterFirstSteps = dataspectrumaftersubstraction;
801 gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
804 printf("after normalize getting IP \n");
806 /////////////////////////////////////////////////////////////////////////////////////////
807 // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
808 /////////////////////////////////////////////////////////////////////////////////////////
810 AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
811 AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
812 AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
814 if(fEfficiencyFunction){
815 dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
816 dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
818 else if(dataContainerV0){
819 dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
820 dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
828 TList *listunfolded = Unfold(dataGridAfterFirstSteps);
830 printf("Unfolded failed\n");
833 THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
834 THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
835 if(!correctedspectrum){
836 AliError("No corrected spectrum\n");
839 if(!residualspectrum){
840 AliError("No residual spectrum\n");
844 /////////////////////
848 AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
855 if(fDebugLevel > 0.0) {
858 if(fBeamType==0) ptpr=0;
859 if(fBeamType==1) ptpr=1;
861 TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
862 ccorrected->Divide(2,1);
865 TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
866 correctedspectrumD->SetTitle("");
867 correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
868 correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
869 correctedspectrumD->SetMarkerStyle(26);
870 correctedspectrumD->SetMarkerColor(kBlue);
871 correctedspectrumD->SetLineColor(kBlue);
872 correctedspectrumD->Draw("AP");
873 TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
874 alltogetherspectrumD->SetTitle("");
875 alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
876 alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
877 alltogetherspectrumD->SetMarkerStyle(25);
878 alltogetherspectrumD->SetMarkerColor(kBlack);
879 alltogetherspectrumD->SetLineColor(kBlack);
880 alltogetherspectrumD->Draw("P");
881 TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
882 legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
883 legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
884 legcorrected->Draw("same");
886 TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
887 TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
888 TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
889 ratiocorrected->SetName("ratiocorrected");
890 ratiocorrected->SetTitle("");
891 ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
892 ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
893 ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
894 ratiocorrected->SetStats(0);
895 ratiocorrected->Draw();
896 if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
900 CalculateNonHFEsyst(0);
904 // Dump to file if needed
907 // to do centrality dependent
910 out = new TFile("finalSpectrum.root","recreate");
913 correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
914 correctedspectrumD->Write();
915 alltogetherspectrumD->SetName("AlltogetherSpectrum");
916 alltogetherspectrumD->Write();
917 ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
918 ratiocorrected->Write();
920 correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
921 correctedspectrum->Write();
922 alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
923 alltogetherCorrection->Write();
925 if(unnormalizedRawSpectrum) {
926 unnormalizedRawSpectrum->SetName("beautyAfterIP");
927 unnormalizedRawSpectrum->Write();
930 if(gNormalizedRawSpectrum){
931 gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
932 gNormalizedRawSpectrum->Write();
937 fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
938 fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
939 fEfficiencyCharmSigD[countpp]->Write();
940 fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
941 fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
942 fEfficiencyBeautySigD[countpp]->Write();
943 fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
944 fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
945 fCharmEff[countpp]->Write();
946 fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
947 fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
948 fBeautyEff[countpp]->Write();
949 fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
950 fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
951 fConversionEff[countpp]->Write();
952 fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
953 fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
954 fNonHFEEff[countpp]->Write();
959 TGraphErrors* correctedspectrumDc[kCentrality];
960 TGraphErrors* alltogetherspectrumDc[kCentrality];
961 for(Int_t i=0;i<kCentrality-2;i++)
963 correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
964 correctedspectrumDc[i] = Normalize(correctedspectrum,i);
965 if(correctedspectrumDc[i]){
966 correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
967 correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
968 correctedspectrumDc[i]->Write();
970 alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
971 alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
972 if(alltogetherspectrumDc[i]){
973 alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
974 alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
975 alltogetherspectrumDc[i]->Write();
978 TH1D *centrcrosscheck = correctedspectrum->Projection(0);
979 centrcrosscheck->SetTitle(Form("centrality_%i",i));
980 centrcrosscheck->SetName(Form("centrality_%i",i));
981 centrcrosscheck->Write();
983 TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
984 TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
986 TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0);
987 centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
988 centrcrosscheck2->SetName(Form("centrality2_%i",i));
989 centrcrosscheck2->Write();
991 TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
992 ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
993 ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
994 ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
995 ratiocorrectedc->Write();
997 fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
998 fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
999 fEfficiencyCharmSigD[i]->Write();
1000 fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
1001 fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
1002 fEfficiencyBeautySigD[i]->Write();
1003 fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
1004 fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
1005 fCharmEff[i]->Write();
1006 fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
1007 fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
1008 fBeautyEff[i]->Write();
1009 fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
1010 fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
1011 fConversionEff[i]->Write();
1012 fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
1013 fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
1014 fNonHFEEff[i]->Write();
1027 //____________________________________________________________
1028 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
1030 // Apply background subtraction
1047 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1049 AliError("Data Container not available");
1052 printf("Step data: %d\n",fStepData);
1053 AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
1055 AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
1056 dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
1059 // Background Estimate
1060 AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
1061 if(!backgroundContainer){
1062 AliError("MC background container not found");
1066 Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
1067 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
1069 if(!fInclusiveSpectrum){
1070 //Background subtraction for IP analysis
1072 TH1D *incElecCent[kCentrality-1];
1073 TH1D *charmCent[kCentrality-1];
1074 TH1D *convCent[kCentrality-1];
1075 TH1D *nonHFECent[kCentrality-1];
1076 TH1D *subtractedCent[kCentrality-1];
1077 TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
1078 CorrectFromTheWidth(measuredTH1Draw);
1080 THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1081 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1082 sparseIncElec->GetAxis(0)->SetRange(icent,icent);
1083 incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
1084 CorrectFromTheWidth(incElecCent[icent-1]);
1087 TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
1089 rawspectra->SetLogy();
1090 gStyle->SetOptStat(0);
1091 TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
1092 measuredTH1Draw->SetMarkerStyle(20);
1093 measuredTH1Draw->Draw();
1094 measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
1095 lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
1097 Int_t* bins=new Int_t[2];
1098 if(fIPanaHadronBgSubtract){
1099 // Hadron background
1100 printf("Hadron background for IP analysis subtracted!\n");
1104 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1105 bins[0]=htemp->GetNbinsX();
1109 htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
1110 bins[0]=htemp->GetNbinsX();
1111 htemp = (TH1D *) fHadronEffbyIPcut->Projection(1);
1112 bins[1]=htemp->GetNbinsX();
1114 AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
1115 hbgContainer->SetGrid(fHadronEffbyIPcut);
1116 backgroundGrid->Multiply(hbgContainer,1);
1117 // draw raw hadron bg spectra
1118 TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
1119 CorrectFromTheWidth(hadronbg);
1120 hadronbg->SetMarkerColor(7);
1121 hadronbg->SetMarkerStyle(20);
1123 hadronbg->Draw("samep");
1124 lRaw->AddEntry(hadronbg,"hadrons");
1125 // subtract hadron contamination
1126 spectrumSubtracted->Add(backgroundGrid,-1.0);
1128 if(fIPanaCharmBgSubtract){
1130 printf("Charm background for IP analysis subtracted!\n");
1131 AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
1132 // draw charm bg spectra
1133 TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
1134 CorrectFromTheWidth(charmbg);
1135 charmbg->SetMarkerColor(3);
1136 charmbg->SetMarkerStyle(20);
1138 charmbg->Draw("samep");
1139 lRaw->AddEntry(charmbg,"charm elecs");
1140 // subtract charm background
1141 spectrumSubtracted->Add(charmbgContainer,-1.0);
1143 THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
1144 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1145 sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
1146 charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
1147 CorrectFromTheWidth(charmCent[icent-1]);
1151 if(fIPanaConversionBgSubtract){
1152 // Conversion background
1153 AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
1154 // draw conversion bg spectra
1155 TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
1156 CorrectFromTheWidth(conversionbg);
1157 conversionbg->SetMarkerColor(4);
1158 conversionbg->SetMarkerStyle(20);
1160 conversionbg->Draw("samep");
1161 lRaw->AddEntry(conversionbg,"conversion elecs");
1162 // subtract conversion background
1163 spectrumSubtracted->Add(conversionbgContainer,-1.0);
1165 THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
1166 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1167 sparseconvElec->GetAxis(0)->SetRange(icent,icent);
1168 convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
1169 CorrectFromTheWidth(convCent[icent-1]);
1173 if(fIPanaNonHFEBgSubtract){
1174 // NonHFE background
1175 AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
1176 // draw Dalitz/dielectron bg spectra
1177 TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
1178 CorrectFromTheWidth(nonhfebg);
1179 nonhfebg->SetMarkerColor(6);
1180 nonhfebg->SetMarkerStyle(20);
1182 nonhfebg->Draw("samep");
1183 lRaw->AddEntry(nonhfebg,"non-HF elecs");
1184 // subtract Dalitz/dielectron background
1185 spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
1187 THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
1188 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1189 sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
1190 nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
1191 CorrectFromTheWidth(nonHFECent[icent-1]);
1196 TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
1197 CorrectFromTheWidth(rawbgsubtracted);
1198 rawbgsubtracted->SetMarkerStyle(24);
1200 lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
1201 rawbgsubtracted->Draw("samep");
1204 //rawspectra->SaveAs("rawspectra.eps");
1207 THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1208 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1209 sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
1210 subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
1211 CorrectFromTheWidth(subtractedCent[icent-1]);
1214 TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
1215 TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
1216 centRaw->Divide(3,3);
1217 for(Int_t icent = 1; icent < kCentrality-1; icent++){
1221 incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
1222 incElecCent[icent-1]->Draw("p");
1223 incElecCent[icent-1]->SetMarkerColor(1);
1224 incElecCent[icent-1]->SetMarkerStyle(20);
1225 charmCent[icent-1]->Draw("samep");
1226 charmCent[icent-1]->SetMarkerColor(3);
1227 charmCent[icent-1]->SetMarkerStyle(20);
1228 convCent[icent-1]->Draw("samep");
1229 convCent[icent-1]->SetMarkerColor(4);
1230 convCent[icent-1]->SetMarkerStyle(20);
1231 nonHFECent[icent-1]->Draw("samep");
1232 nonHFECent[icent-1]->SetMarkerColor(6);
1233 nonHFECent[icent-1]->SetMarkerStyle(20);
1234 subtractedCent[icent-1]->Draw("samep");
1235 subtractedCent[icent-1]->SetMarkerStyle(24);
1237 lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
1238 lCentRaw->AddEntry(charmCent[0],"charm elecs");
1239 lCentRaw->AddEntry(convCent[0],"conversion elecs");
1240 lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
1241 lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
1242 lCentRaw->Draw("SAME");
1252 spectrumSubtracted->Add(backgroundGrid,-1.0);
1256 if(fBackground) delete fBackground;
1257 fBackground = backgroundGrid;
1258 } else delete backgroundGrid;
1261 if(fDebugLevel > 0) {
1264 if(fBeamType==0) ptprd=0;
1265 if(fBeamType==1) ptprd=1;
1267 TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
1268 cbackgroundsubtraction->Divide(3,1);
1269 cbackgroundsubtraction->cd(1);
1271 TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
1272 TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
1273 CorrectFromTheWidth(measuredTH1Daftersubstraction);
1274 CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
1275 measuredTH1Daftersubstraction->SetStats(0);
1276 measuredTH1Daftersubstraction->SetTitle("");
1277 measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1278 measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1279 measuredTH1Daftersubstraction->SetMarkerStyle(25);
1280 measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
1281 measuredTH1Daftersubstraction->SetLineColor(kBlack);
1282 measuredTH1Dbeforesubstraction->SetStats(0);
1283 measuredTH1Dbeforesubstraction->SetTitle("");
1284 measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1285 measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1286 measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
1287 measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
1288 measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
1289 measuredTH1Daftersubstraction->Draw();
1290 measuredTH1Dbeforesubstraction->Draw("same");
1291 TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
1292 legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
1293 legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
1294 legsubstraction->Draw("same");
1295 cbackgroundsubtraction->cd(2);
1297 TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1298 ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1299 ratiomeasuredcontamination->SetTitle("");
1300 ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1301 ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1302 ratiomeasuredcontamination->Sumw2();
1303 ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1304 ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1305 ratiomeasuredcontamination->SetStats(0);
1306 ratiomeasuredcontamination->SetMarkerStyle(26);
1307 ratiomeasuredcontamination->SetMarkerColor(kBlack);
1308 ratiomeasuredcontamination->SetLineColor(kBlack);
1309 for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1310 ratiomeasuredcontamination->SetBinError(k+1,0.0);
1312 ratiomeasuredcontamination->Draw("P");
1313 cbackgroundsubtraction->cd(3);
1314 TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
1315 CorrectFromTheWidth(measuredTH1background);
1316 measuredTH1background->SetStats(0);
1317 measuredTH1background->SetTitle("");
1318 measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1319 measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1320 measuredTH1background->SetMarkerStyle(26);
1321 measuredTH1background->SetMarkerColor(kBlack);
1322 measuredTH1background->SetLineColor(kBlack);
1323 measuredTH1background->Draw();
1324 if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
1328 TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1329 cbackgrounde->Divide(2,1);
1330 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1331 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1333 THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1334 TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1335 THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1336 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1337 Int_t nbbin = cenaxisb->GetNbins();
1338 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1339 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1340 for(Int_t binc = 0; binc < nbbin; binc++){
1341 TString titlee("BackgroundSubtraction_centrality_bin_");
1343 TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1344 cbackground->Divide(2,1);
1347 cenaxisa->SetRange(binc+1,binc+1);
1348 cenaxisb->SetRange(binc+1,binc+1);
1349 TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1350 TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1351 CorrectFromTheWidth(aftersubstraction);
1352 CorrectFromTheWidth(beforesubstraction);
1353 aftersubstraction->SetStats(0);
1354 aftersubstraction->SetTitle((const char*)titlee);
1355 aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1356 aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1357 aftersubstraction->SetMarkerStyle(25);
1358 aftersubstraction->SetMarkerColor(kBlack);
1359 aftersubstraction->SetLineColor(kBlack);
1360 beforesubstraction->SetStats(0);
1361 beforesubstraction->SetTitle((const char*)titlee);
1362 beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1363 beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1364 beforesubstraction->SetMarkerStyle(24);
1365 beforesubstraction->SetMarkerColor(kBlue);
1366 beforesubstraction->SetLineColor(kBlue);
1367 aftersubstraction->Draw();
1368 beforesubstraction->Draw("same");
1369 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1370 lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1371 lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1373 cbackgrounde->cd(1);
1375 TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1376 aftersubtractionn->SetMarkerStyle(stylee[binc]);
1377 aftersubtractionn->SetMarkerColor(colorr[binc]);
1378 if(binc==0) aftersubtractionn->Draw();
1379 else aftersubtractionn->Draw("same");
1380 legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1381 cbackgrounde->cd(2);
1383 TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1384 aftersubtractionng->SetMarkerStyle(stylee[binc]);
1385 aftersubtractionng->SetMarkerColor(colorr[binc]);
1386 if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1387 if(binc==0) aftersubtractionng->Draw();
1388 else aftersubtractionng->Draw("same");
1389 legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1391 TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1392 ratiocontamination->SetName("ratiocontamination");
1393 ratiocontamination->SetTitle((const char*)titlee);
1394 ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1395 ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1396 ratiocontamination->Add(aftersubstraction,-1.0);
1397 ratiocontamination->Divide(beforesubstraction);
1398 Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1399 for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1400 ratiocontamination->SetBinError(nbinpt+1,0.0);
1402 ratiocontamination->SetStats(0);
1403 ratiocontamination->SetMarkerStyle(26);
1404 ratiocontamination->SetMarkerColor(kBlack);
1405 ratiocontamination->SetLineColor(kBlack);
1406 ratiocontamination->Draw("P");
1409 cbackgrounde->cd(1);
1410 legtotal->Draw("same");
1411 cbackgrounde->cd(2);
1412 legtotalg->Draw("same");
1414 cenaxisa->SetRange(0,nbbin);
1415 cenaxisb->SetRange(0,nbbin);
1416 if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
1420 return spectrumSubtracted;
1423 //____________________________________________________________
1424 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1426 // calculate charm background
1441 if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
1443 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1445 AliError("MC Container not available");
1450 AliError("No Correlation map available");
1454 AliCFDataGrid *charmBackgroundGrid= 0x0;
1455 charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1457 TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1458 Int_t* bins=new Int_t[2];
1459 bins[0]=charmbgaftertofpid->GetNbinsX();
1464 charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
1465 bins[1]=charmbgaftertofpid->GetNbinsX();
1469 AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
1470 ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1471 TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
1473 charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
1475 Int_t *nBinpp=new Int_t[1];
1476 Int_t* binspp=new Int_t[1];
1477 binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
1479 Int_t *nBinPbPb=new Int_t[2];
1480 Int_t* binsPbPb=new Int_t[2];
1481 binsPbPb[1]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
1484 Int_t looppt=binspp[0];
1485 if(fBeamType==1) looppt=binsPbPb[1];
1487 for(Long_t iBin=1; iBin<= looppt;iBin++){
1491 charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
1492 charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
1496 // loop over centrality
1497 for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1500 Double_t evtnormPbPb=0;
1501 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1502 charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1503 charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1508 TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
1510 AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
1511 weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1512 TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
1513 charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1514 TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
1516 // Efficiency (set efficiency to 1 for only folding)
1517 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1518 efficiencyD->CalculateEfficiency(0,0);
1521 if(fBeamType==0)nDim = 1;
1522 if(fBeamType==1)nDim = 2;
1523 AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1524 folding.SetMaxNumberOfIterations(1);
1528 THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1529 THnSparse* result=(THnSparse*)result1->Clone();
1530 charmBackgroundGrid->SetGrid(result);
1531 TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
1533 //Charm background evaluation plots
1535 TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1536 cCharmBgEval->Divide(3,1);
1538 cCharmBgEval->cd(1);
1540 if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
1543 Double_t evtnormPbPb=0;
1544 for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1546 if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1548 charmbgaftertofpid->Scale(evtnormPbPb);
1551 CorrectFromTheWidth(charmbgaftertofpid);
1552 charmbgaftertofpid->SetMarkerStyle(25);
1553 charmbgaftertofpid->Draw("p");
1554 charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1555 charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1558 CorrectFromTheWidth(charmbgafteripcut);
1559 charmbgafteripcut->SetMarkerStyle(24);
1560 charmbgafteripcut->Draw("samep");
1562 CorrectFromTheWidth(charmbgafterweight);
1563 charmbgafterweight->SetMarkerStyle(24);
1564 charmbgafterweight->SetMarkerColor(4);
1565 charmbgafterweight->Draw("samep");
1567 CorrectFromTheWidth(charmbgafterfolding);
1568 charmbgafterfolding->SetMarkerStyle(24);
1569 charmbgafterfolding->SetMarkerColor(2);
1570 charmbgafterfolding->Draw("samep");
1572 cCharmBgEval->cd(2);
1573 parametrizedcharmpidipeff->SetMarkerStyle(24);
1574 parametrizedcharmpidipeff->Draw("p");
1575 parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1577 cCharmBgEval->cd(3);
1578 charmweightingfc->SetMarkerStyle(24);
1579 charmweightingfc->Draw("p");
1580 charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1581 charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1583 cCharmBgEval->cd(1);
1584 TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1585 legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1586 legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1587 legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1588 legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1589 legcharmbg->Draw("same");
1591 cCharmBgEval->cd(2);
1592 TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1593 legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1594 legcharmbg2->Draw("same");
1596 CorrectStatErr(charmBackgroundGrid);
1597 if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
1605 if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
1607 return charmBackgroundGrid;
1610 //____________________________________________________________
1611 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1613 // calculate conversion background
1616 Double_t evtnorm[1] = {0.0};
1617 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1618 printf("check event!!! %lf \n",evtnorm[0]);
1620 AliCFContainer *backgroundContainer = 0x0;
1623 backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
1624 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1625 backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
1629 // Background Estimate
1630 backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1632 if(!backgroundContainer){
1633 AliError("MC background container not found");
1637 Int_t stepbackground = 3;
1639 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1640 Int_t *nBinpp=new Int_t[1];
1641 Int_t* binspp=new Int_t[1];
1642 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1644 Int_t *nBinPbPb=new Int_t[2];
1645 Int_t* binsPbPb=new Int_t[2];
1646 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1649 Int_t looppt=binspp[0];
1650 if(fBeamType==1) looppt=binsPbPb[1];
1652 for(Long_t iBin=1; iBin<= looppt;iBin++){
1656 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1657 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1661 // loop over centrality
1662 for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1665 Double_t evtnormPbPb=0;
1666 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1667 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1668 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1672 //end of workaround for statistical errors
1674 AliCFDataGrid *weightedConversionContainer;
1675 if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
1676 else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
1677 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1678 backgroundGrid->Multiply(weightedConversionContainer,1.0);
1685 return backgroundGrid;
1689 //____________________________________________________________
1690 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1692 // calculate non-HFE background
1695 Double_t evtnorm[1] = {0.0};
1696 if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1697 printf("check event!!! %lf \n",evtnorm[0]);
1699 AliCFContainer *backgroundContainer = 0x0;
1701 backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
1702 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1704 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
1706 backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
1710 // Background Estimate
1711 backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1713 if(!backgroundContainer){
1714 AliError("MC background container not found");
1719 Int_t stepbackground = 3;
1721 AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1722 Int_t *nBinpp=new Int_t[1];
1723 Int_t* binspp=new Int_t[1];
1724 binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1726 Int_t *nBinPbPb=new Int_t[2];
1727 Int_t* binsPbPb=new Int_t[2];
1728 binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
1731 Int_t looppt=binspp[0];
1732 if(fBeamType==1) looppt=binsPbPb[1];
1735 for(Long_t iBin=1; iBin<= looppt;iBin++){
1739 backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1740 backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1744 for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1747 Double_t evtnormPbPb=0;
1748 if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1749 backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1750 backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1754 //end of workaround for statistical errors
1755 AliCFDataGrid *weightedNonHFEContainer;
1756 if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1757 else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
1758 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1759 backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1766 return backgroundGrid;
1769 //____________________________________________________________
1770 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1773 // Apply TPC pid efficiency correction from parametrisation
1776 // Data in the right format
1777 AliCFDataGrid *dataGrid = 0x0;
1779 dataGrid = bgsubpectrum;
1783 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1785 AliError("Data Container not available");
1788 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1790 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1791 result->SetName("ParametrizedEfficiencyBefore");
1792 THnSparse *h = result->GetGrid();
1793 Int_t nbdimensions = h->GetNdimensions();
1794 //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1795 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1797 AliError("Data Container not available");
1800 AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1801 dataContainerbis->Add(dataContainerbis,-1.0);
1804 Int_t* coord = new Int_t[nbdimensions];
1805 memset(coord, 0, sizeof(Int_t) * nbdimensions);
1806 Double_t* points = new Double_t[nbdimensions];
1808 ULong64_t nEntries = h->GetNbins();
1809 for (ULong64_t i = 0; i < nEntries; ++i) {
1811 Double_t value = h->GetBinContent(i, coord);
1812 //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1813 //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1815 // Get the bin co-ordinates given an coord
1816 for (Int_t j = 0; j < nbdimensions; ++j)
1817 points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1819 if (!fEfficiencyFunction->IsInside(points))
1821 TF1::RejectPoint(kFALSE);
1823 // Evaulate function at points
1824 Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1825 //printf("Value efficiency is %f\n",valueEfficiency);
1827 if(valueEfficiency > 0.0) {
1828 h->SetBinContent(coord,value/valueEfficiency);
1829 dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1831 Double_t error = h->GetBinError(i);
1832 h->SetBinError(coord,error/valueEfficiency);
1833 dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1841 AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1843 if(fDebugLevel > 0) {
1845 TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1846 cEfficiencyParametrized->Divide(2,1);
1847 cEfficiencyParametrized->cd(1);
1848 TH1D *afterE = (TH1D *) resultt->Project(0);
1849 TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1850 CorrectFromTheWidth(afterE);
1851 CorrectFromTheWidth(beforeE);
1852 afterE->SetStats(0);
1853 afterE->SetTitle("");
1854 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1855 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1856 afterE->SetMarkerStyle(25);
1857 afterE->SetMarkerColor(kBlack);
1858 afterE->SetLineColor(kBlack);
1859 beforeE->SetStats(0);
1860 beforeE->SetTitle("");
1861 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1862 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1863 beforeE->SetMarkerStyle(24);
1864 beforeE->SetMarkerColor(kBlue);
1865 beforeE->SetLineColor(kBlue);
1868 beforeE->Draw("same");
1869 TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1870 legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1871 legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1872 legefficiencyparametrized->Draw("same");
1873 cEfficiencyParametrized->cd(2);
1874 fEfficiencyFunction->Draw();
1875 //cEfficiencyParametrized->cd(3);
1876 //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1877 //ratioefficiency->Divide(afterE);
1878 //ratioefficiency->Draw();
1880 if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
1887 //____________________________________________________________
1888 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1891 // Apply TPC pid efficiency correction from V0
1894 AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1896 AliError("V0 Container not available");
1901 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1902 efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1904 // Data in the right format
1905 AliCFDataGrid *dataGrid = 0x0;
1907 dataGrid = bgsubpectrum;
1911 AliCFContainer *dataContainer = GetContainer(kDataContainer);
1913 AliError("Data Container not available");
1917 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1921 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1922 result->ApplyEffCorrection(*efficiencyD);
1924 if(fDebugLevel > 0) {
1927 if(fBeamType==0) ptpr=0;
1928 if(fBeamType==1) ptpr=1;
1930 TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1931 cV0Efficiency->Divide(2,1);
1932 cV0Efficiency->cd(1);
1933 TH1D *afterE = (TH1D *) result->Project(ptpr);
1934 TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1935 afterE->SetStats(0);
1936 afterE->SetTitle("");
1937 afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1938 afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1939 afterE->SetMarkerStyle(25);
1940 afterE->SetMarkerColor(kBlack);
1941 afterE->SetLineColor(kBlack);
1942 beforeE->SetStats(0);
1943 beforeE->SetTitle("");
1944 beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1945 beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1946 beforeE->SetMarkerStyle(24);
1947 beforeE->SetMarkerColor(kBlue);
1948 beforeE->SetLineColor(kBlue);
1950 beforeE->Draw("same");
1951 TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1952 legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1953 legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1954 legV0efficiency->Draw("same");
1955 cV0Efficiency->cd(2);
1956 TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1957 efficiencyDproj->SetTitle("");
1958 efficiencyDproj->SetStats(0);
1959 efficiencyDproj->SetMarkerStyle(25);
1960 efficiencyDproj->Draw();
1964 TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1965 cV0Efficiencye->Divide(2,1);
1966 TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1967 TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1969 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1970 TAxis *cenaxisa = sparseafter->GetAxis(0);
1971 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1972 TAxis *cenaxisb = sparsebefore->GetAxis(0);
1973 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1974 TAxis *cenaxisc = efficiencya->GetAxis(0);
1975 Int_t nbbin = cenaxisb->GetNbins();
1976 Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1977 Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1978 for(Int_t binc = 0; binc < nbbin; binc++){
1979 TString titlee("V0Efficiency_centrality_bin_");
1981 TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1982 ccV0Efficiency->Divide(2,1);
1983 ccV0Efficiency->cd(1);
1985 cenaxisa->SetRange(binc+1,binc+1);
1986 cenaxisb->SetRange(binc+1,binc+1);
1987 cenaxisc->SetRange(binc+1,binc+1);
1988 TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1989 TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1990 CorrectFromTheWidth(aftere);
1991 CorrectFromTheWidth(beforee);
1992 aftere->SetStats(0);
1993 aftere->SetTitle((const char*)titlee);
1994 aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1995 aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1996 aftere->SetMarkerStyle(25);
1997 aftere->SetMarkerColor(kBlack);
1998 aftere->SetLineColor(kBlack);
1999 beforee->SetStats(0);
2000 beforee->SetTitle((const char*)titlee);
2001 beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2002 beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2003 beforee->SetMarkerStyle(24);
2004 beforee->SetMarkerColor(kBlue);
2005 beforee->SetLineColor(kBlue);
2007 beforee->Draw("same");
2008 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2009 lega->AddEntry(beforee,"Before correction","p");
2010 lega->AddEntry(aftere,"After correction","p");
2012 cV0Efficiencye->cd(1);
2014 TH1D *afteree = (TH1D *) aftere->Clone();
2015 afteree->SetMarkerStyle(stylee[binc]);
2016 afteree->SetMarkerColor(colorr[binc]);
2017 if(binc==0) afteree->Draw();
2018 else afteree->Draw("same");
2019 legtotal->AddEntry(afteree,(const char*) titlee,"p");
2020 cV0Efficiencye->cd(2);
2022 TH1D *aftereeu = (TH1D *) aftere->Clone();
2023 aftereeu->SetMarkerStyle(stylee[binc]);
2024 aftereeu->SetMarkerColor(colorr[binc]);
2025 if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
2026 if(binc==0) aftereeu->Draw();
2027 else aftereeu->Draw("same");
2028 legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
2029 ccV0Efficiency->cd(2);
2030 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
2031 efficiencyDDproj->SetTitle("");
2032 efficiencyDDproj->SetStats(0);
2033 efficiencyDDproj->SetMarkerStyle(25);
2034 efficiencyDDproj->Draw();
2037 cV0Efficiencye->cd(1);
2038 legtotal->Draw("same");
2039 cV0Efficiencye->cd(2);
2040 legtotalg->Draw("same");
2042 cenaxisa->SetRange(0,nbbin);
2043 cenaxisb->SetRange(0,nbbin);
2044 cenaxisc->SetRange(0,nbbin);
2046 if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
2055 //____________________________________________________________
2056 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
2059 // Unfold and eventually correct for efficiency the bgsubspectrum
2062 AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2064 AliError("MC Container not available");
2069 AliError("No Correlation map available");
2074 AliCFDataGrid *dataGrid = 0x0;
2076 dataGrid = bgsubpectrum;
2080 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2082 AliError("Data Container not available");
2086 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2090 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2091 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2094 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2095 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2097 if(!fBeauty2ndMethod)
2099 // Consider parameterized IP cut efficiency
2100 if(!fInclusiveSpectrum){
2101 Int_t* bins=new Int_t[1];
2104 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2105 beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
2106 efficiencyD->Multiply(beffContainer,1);
2113 AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2114 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2115 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2116 if(fNRandomIter > 0) unfolding.SetNRandomIterations(fNRandomIter);
2117 if(fSetSmoothing) unfolding.UseSmoothing();
2121 THnSparse* result = unfolding.GetUnfolded();
2122 THnSparse* residual = unfolding.GetEstMeasured();
2123 TList *listofresults = new TList;
2124 listofresults->SetOwner();
2125 listofresults->AddAt((THnSparse*)result->Clone(),0);
2126 listofresults->AddAt((THnSparse*)residual->Clone(),1);
2128 if(fDebugLevel > 0) {
2131 if(fBeamType==0) ptpr=0;
2132 if(fBeamType==1) ptpr=1;
2134 TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
2135 cresidual->Divide(2,1);
2138 TGraphErrors* residualspectrumD = Normalize(residual);
2139 if(!residualspectrumD) {
2140 AliError("Number of Events not set for the normalization");
2143 residualspectrumD->SetTitle("");
2144 residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2145 residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2146 residualspectrumD->SetMarkerStyle(26);
2147 residualspectrumD->SetMarkerColor(kBlue);
2148 residualspectrumD->SetLineColor(kBlue);
2149 residualspectrumD->Draw("AP");
2150 AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2151 dataGridBis->SetName("dataGridBis");
2152 TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2153 measuredspectrumD->SetTitle("");
2154 measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2155 measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2156 measuredspectrumD->SetMarkerStyle(25);
2157 measuredspectrumD->SetMarkerColor(kBlack);
2158 measuredspectrumD->SetLineColor(kBlack);
2159 measuredspectrumD->Draw("P");
2160 TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2161 legres->AddEntry(residualspectrumD,"Residual","p");
2162 legres->AddEntry(measuredspectrumD,"Measured","p");
2163 legres->Draw("same");
2165 TH1D *residualTH1D = residual->Projection(ptpr);
2166 TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
2167 TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2168 ratioresidual->SetName("ratioresidual");
2169 ratioresidual->SetTitle("");
2170 ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2171 ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2172 ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2173 ratioresidual->SetStats(0);
2174 ratioresidual->Draw();
2176 if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
2179 return listofresults;
2183 //____________________________________________________________
2184 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
2187 // Apply unfolding and efficiency correction together to bgsubspectrum
2190 AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
2192 AliError("MC Container not available");
2197 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2198 efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2201 if(!fBeauty2ndMethod)
2203 // Consider parameterized IP cut efficiency
2204 if(!fInclusiveSpectrum){
2205 Int_t* bins=new Int_t[1];
2208 AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2209 beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
2210 efficiencyD->Multiply(beffContainer,1);
2214 // Data in the right format
2215 AliCFDataGrid *dataGrid = 0x0;
2217 dataGrid = bgsubpectrum;
2221 AliCFContainer *dataContainer = GetContainer(kDataContainer);
2223 AliError("Data Container not available");
2227 dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2231 AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2232 result->ApplyEffCorrection(*efficiencyD);
2234 if(fDebugLevel > 0) {
2237 if(fBeamType==0) ptpr=0;
2238 if(fBeamType==1) ptpr=1;
2240 printf("Step MC: %d\n",fStepMC);
2241 printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2242 printf("Step MC true: %d\n",fStepTrue);
2243 AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2244 AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2245 AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2247 AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2249 TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2251 TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
2252 efficiencymcPIDD->SetTitle("");
2253 efficiencymcPIDD->SetStats(0);
2254 efficiencymcPIDD->SetMarkerStyle(25);
2255 efficiencymcPIDD->Draw();
2256 TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
2257 efficiencymctrackinggeoD->SetTitle("");
2258 efficiencymctrackinggeoD->SetStats(0);
2259 efficiencymctrackinggeoD->SetMarkerStyle(26);
2260 efficiencymctrackinggeoD->Draw("same");
2261 TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
2262 efficiencymcallD->SetTitle("");
2263 efficiencymcallD->SetStats(0);
2264 efficiencymcallD->SetMarkerStyle(27);
2265 efficiencymcallD->Draw("same");
2266 TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
2267 efficiencyesdallD->SetTitle("");
2268 efficiencyesdallD->SetStats(0);
2269 efficiencyesdallD->SetMarkerStyle(24);
2270 efficiencyesdallD->Draw("same");
2271 TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2272 legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2273 legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2274 legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2275 legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2276 legeff->Draw("same");
2280 THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2281 TAxis *cenaxisa = sparseafter->GetAxis(0);
2282 THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2283 TAxis *cenaxisb = sparsebefore->GetAxis(0);
2284 THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2285 TAxis *cenaxisc = efficiencya->GetAxis(0);
2286 //Int_t nbbin = cenaxisb->GetNbins();
2287 //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2288 //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2289 for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
2290 TString titlee("Efficiency_centrality_bin_");
2291 titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2293 titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2294 TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2295 cefficiencye->Divide(2,1);
2296 cefficiencye->cd(1);
2298 cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2299 cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2300 TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2301 TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2302 CorrectFromTheWidth(afterefficiency);
2303 CorrectFromTheWidth(beforeefficiency);
2304 afterefficiency->SetStats(0);
2305 afterefficiency->SetTitle((const char*)titlee);
2306 afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2307 afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2308 afterefficiency->SetMarkerStyle(25);
2309 afterefficiency->SetMarkerColor(kBlack);
2310 afterefficiency->SetLineColor(kBlack);
2311 beforeefficiency->SetStats(0);
2312 beforeefficiency->SetTitle((const char*)titlee);
2313 beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2314 beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2315 beforeefficiency->SetMarkerStyle(24);
2316 beforeefficiency->SetMarkerColor(kBlue);
2317 beforeefficiency->SetLineColor(kBlue);
2318 afterefficiency->Draw();
2319 beforeefficiency->Draw("same");
2320 TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2321 lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2322 lega->AddEntry(afterefficiency,"After efficiency correction","p");
2324 cefficiencye->cd(2);
2325 cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2326 TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2327 efficiencyDDproj->SetTitle("");
2328 efficiencyDDproj->SetStats(0);
2329 efficiencyDDproj->SetMarkerStyle(25);
2330 efficiencyDDproj->SetMarkerColor(2);
2331 efficiencyDDproj->Draw();
2332 cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2333 TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2334 efficiencyDDproja->SetTitle("");
2335 efficiencyDDproja->SetStats(0);
2336 efficiencyDDproja->SetMarkerStyle(26);
2337 efficiencyDDproja->SetMarkerColor(4);
2338 efficiencyDDproja->Draw("same");
2342 if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
2349 //__________________________________________________________________________________
2350 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
2352 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2353 // Give the final pt spectrum to be compared
2356 if(fNEvents[i] > 0) {
2359 if(fBeamType==0) ptpr=0;
2360 if(fBeamType==1) ptpr=1;
2362 TH1D* projection = spectrum->Projection(ptpr);
2363 CorrectFromTheWidth(projection);
2364 TGraphErrors *graphError = NormalizeTH1(projection,i);
2373 //__________________________________________________________________________________
2374 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
2376 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2377 // Give the final pt spectrum to be compared
2379 if(fNEvents[i] > 0) {
2382 if(fBeamType==0) ptpr=0;
2383 if(fBeamType==1) ptpr=1;
2385 TH1D* projection = (TH1D *) spectrum->Project(ptpr);
2386 CorrectFromTheWidth(projection);
2387 TGraphErrors *graphError = NormalizeTH1(projection,i);
2397 //__________________________________________________________________________________
2398 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2400 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2401 // Give the final pt spectrum to be compared
2403 Double_t chargecoefficient = 0.5;
2404 if(fChargeChoosen != 0) chargecoefficient = 1.0;
2406 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2407 printf("Normalizing Eta Range %f\n", etarange);
2408 if(fNEvents[i] > 0) {
2410 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2411 Double_t p = 0, dp = 0; Int_t point = 1;
2412 Double_t n = 0, dN = 0;
2413 Double_t nCorr = 0, dNcorr = 0;
2414 //Double_t errdN = 0, errdp = 0;
2416 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2417 point = ibin - input->GetXaxis()->GetFirst();
2418 p = input->GetXaxis()->GetBinCenter(ibin);
2419 dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2420 n = input->GetBinContent(ibin);
2421 dN = input->GetBinError(ibin);
2423 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
2424 errdN = 1./(2. * TMath::Pi() * p);
2425 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2426 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2427 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN);
2429 spectrumNormalized->SetPoint(point, p, nCorr);
2430 spectrumNormalized->SetPointError(point, dp, dNcorr);
2432 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2433 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2434 spectrumNormalized->SetMarkerStyle(22);
2435 spectrumNormalized->SetMarkerColor(kBlue);
2436 spectrumNormalized->SetLineColor(kBlue);
2438 return spectrumNormalized;
2446 //__________________________________________________________________________________
2447 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
2449 // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2450 // Give the final pt spectrum to be compared
2452 Double_t chargecoefficient = 0.5;
2453 if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
2455 Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2456 printf("Normalizing Eta Range %f\n", etarange);
2457 if(normalization > 0) {
2459 TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2460 Double_t p = 0, dp = 0; Int_t point = 1;
2461 Double_t n = 0, dN = 0;
2462 Double_t nCorr = 0, dNcorr = 0;
2463 //Double_t errdN = 0, errdp = 0;
2465 for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2466 point = ibin - input->GetXaxis()->GetFirst();
2467 p = input->GetXaxis()->GetBinCenter(ibin);
2468 //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2469 n = input->GetBinContent(ibin);
2470 dN = input->GetBinError(ibin);
2472 nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
2473 errdN = 1./(2. * TMath::Pi() * p);
2474 //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2475 //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2476 dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN);
2478 spectrumNormalized->SetPoint(point, p, nCorr);
2479 spectrumNormalized->SetPointError(point, dp, dNcorr);
2481 spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2482 spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2483 spectrumNormalized->SetMarkerStyle(22);
2484 spectrumNormalized->SetMarkerColor(kBlue);
2485 spectrumNormalized->SetLineColor(kBlue);
2487 return spectrumNormalized;
2495 //____________________________________________________________
2496 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2498 // Set the container for a given step to the
2500 if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
2501 fCFContainers->AddAt(cont, type);
2504 //____________________________________________________________
2505 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2507 // Get Correction Framework Container for given type
2509 if(!fCFContainers) return NULL;
2510 return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2512 //____________________________________________________________
2513 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
2515 // Slice bin for a given source of electron
2516 // nDim is the number of dimension the corrections are done
2517 // dimensions are the definition of the dimensions
2518 // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2519 // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2520 // centrality (-1 means we do not cut on centrality)
2523 Double_t *varMin = new Double_t[container->GetNVar()],
2524 *varMax = new Double_t[container->GetNVar()];
2526 Double_t *binLimits;
2527 for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2529 binLimits = new Double_t[container->GetNBins(ivar)+1];
2530 container->GetBinLimits(ivar,binLimits);
2531 varMin[ivar] = binLimits[0];
2532 varMax[ivar] = binLimits[container->GetNBins(ivar)];
2535 if((source>= 0) && (source<container->GetNBins(ivar))) {
2536 varMin[ivar] = container->GetAxis(4,0)->GetBinLowEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2537 varMax[ivar] = container->GetAxis(4,0)->GetBinUpEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2542 if(charge != kAllCharge){
2543 varMin[ivar] = container->GetAxis(3,0)->GetBinLowEdge(container->GetAxis(3,0)->FindBin(charge));
2544 varMax[ivar] = container->GetAxis(3,0)->GetBinUpEdge(container->GetAxis(3,0)->FindBin(charge));
2549 for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2550 AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
2552 varMin[ivar] = fEtaRange[0];
2553 varMax[ivar] = fEtaRange[1];
2557 fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2558 fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2559 AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2562 if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
2563 varMin[ivar] = binLimits[centralitylow];
2564 varMax[ivar] = binLimits[centralityhigh];
2566 TAxis *axistest = container->GetAxis(5,0);
2567 printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
2568 printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
2569 Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
2570 Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
2571 printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
2581 AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2582 AddTemporaryObject(k);
2583 delete[] varMin; delete[] varMax;
2589 //_________________________________________________________________________
2590 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
2592 // Slice correlation
2595 Int_t ndimensions = correlationmatrix->GetNdimensions();
2596 //printf("Number of dimension %d correlation map\n",ndimensions);
2597 if(ndimensions < (2*nDim)) {
2598 AliError("Problem in the dimensions");
2602 // Cut in centrality is centrality > -1
2603 if((centralitylow >=0) && (centralityhigh >=0)) {
2605 TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
2606 TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
2608 Int_t bins0 = axiscentrality0->GetNbins();
2609 Int_t bins1 = axiscentrality1->GetNbins();
2611 printf("Number of centrality bins: %d and %d\n",bins0,bins1);
2612 if(bins0 != bins1) {
2613 AliError("Problem in the dimensions");
2617 if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
2618 axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
2619 axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
2621 Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
2622 Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
2623 Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
2624 Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
2625 printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
2626 printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
2628 TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
2629 ctestcorrelation->cd(1);
2630 TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
2631 projection->Draw("colz");
2638 Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2639 //printf("Number of dimension %d container\n",ndimensionsContainer);
2641 Int_t *dim = new Int_t[nDim*2];
2642 for(Int_t iter=0; iter < nDim; iter++){
2643 dim[iter] = dimensions[iter];
2644 dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2645 //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2648 THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2654 //___________________________________________________________________________
2655 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2657 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2660 TAxis *axis = h1->GetXaxis();
2661 Int_t nbinX = h1->GetNbinsX();
2663 for(Int_t i = 1; i <= nbinX; i++) {
2665 Double_t width = axis->GetBinWidth(i);
2666 Double_t content = h1->GetBinContent(i);
2667 Double_t error = h1->GetBinError(i);
2668 h1->SetBinContent(i,content/width);
2669 h1->SetBinError(i,error/width);
2674 //___________________________________________________________________________
2675 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2677 // Correct statistical error
2680 TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2681 Int_t nbinX = h1->GetNbinsX();
2683 for(Long_t i = 1; i <= nbinX; i++) {
2685 Float_t content = h1->GetBinContent(i);
2687 Float_t error = TMath::Sqrt(content);
2688 backgroundGrid->SetElementError(bins, error);
2693 //____________________________________________________________
2694 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2696 // Emulate garbage collection on class level
2698 if(!fTemporaryObjects) {
2699 fTemporaryObjects= new TList;
2700 fTemporaryObjects->SetOwner();
2702 fTemporaryObjects->Add(o);
2705 //____________________________________________________________
2706 void AliHFEspectrum::ClearObject(TObject *o){
2708 // Do a safe deletion
2710 if(fTemporaryObjects){
2711 if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2718 //_________________________________________________________________________
2719 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2720 AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2723 //_________________________________________________________________________
2724 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2726 // Create efficiency grid and calculate efficiency
2729 TString name("eff");
2732 AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2733 eff->CalculateEfficiency(step,step0);
2737 //____________________________________________________________________________
2738 THnSparse* AliHFEspectrum::GetCharmWeights(){
2741 // Measured D->e based weighting factors
2744 const Int_t nDimpp=1;
2745 Int_t nBinpp[nDimpp] = {35};
2746 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.};
2747 const Int_t nDimPbPb=2;
2748 Int_t nBinPbPb[nDimPbPb] = {11,35};
2749 Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2751 Int_t looppt=nBinpp[0];
2754 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2755 fWeightCharm->SetBinEdges(0, ptbinning1);
2759 fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2760 fWeightCharm->SetBinEdges(1, ptbinning1);
2761 fWeightCharm->SetBinEdges(0, kCentralityRange);
2762 loopcentr=nBinPbPb[0];
2765 // Weighting factor for pp
2766 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};
2768 // Weighting factor for PbPb (0-20%)
2769 //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};
2771 // Weighting factor for PbPb (40-80%)
2772 //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};
2776 Double_t contents[2];
2778 for(int icentr=0; icentr<loopcentr; icentr++)
2780 for(int i=0; i<looppt; i++){
2781 pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2792 fWeightCharm->Fill(contents,weight[i]);
2796 Int_t nDimSparse = fWeightCharm->GetNdimensions();
2797 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2798 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2800 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2801 binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2802 nBins *= binsvar[iVar];
2805 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2806 // loop that sets 0 error in each bin
2807 for (Long_t iBin=0; iBin<nBins; iBin++) {
2808 Long_t bintmp = iBin ;
2809 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2810 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2811 bintmp /= binsvar[iVar] ;
2813 fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
2819 return fWeightCharm;
2822 //____________________________________________________________________________
2823 void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
2825 // TOF PID efficiencies
2827 if(fBeamType==0) ptpr=0;
2828 if(fBeamType==1) ptpr=1;
2831 const Int_t nCentralitybinning=11; //number of centrality bins
2834 loopcentr=nCentralitybinning;
2837 TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2838 TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2839 TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
2841 TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
2842 cefficiencyParamtof->cd();
2844 AliCFContainer *mccontainermcD = 0x0;
2845 AliCFContainer *mccontaineresdD = 0x0;
2846 TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2847 TH1D* efficiencyTOFPIDD[nCentralitybinning];
2848 TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2849 TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
2850 Int_t source = -1; //get parameterized TOF PID efficiencies
2852 for(int icentr=0; icentr<loopcentr; icentr++) {
2854 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2855 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2856 AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2857 efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2859 // mb sample for double check
2860 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2861 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2862 AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2863 efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2865 // mb sample with reconstructed variables
2866 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2867 else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2868 AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2869 efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2871 // mb sample with reconstructed variables
2872 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2873 else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2874 AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2875 efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2878 efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2879 efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
2880 efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2881 efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
2882 efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2883 efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
2884 efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2885 efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
2887 //fit (mc enhenced sample)
2888 fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2889 efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2890 efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2891 fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
2893 //fit (esd enhenced sample)
2894 efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2895 efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2896 fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2900 // draw (for PbPb, only 1st bin)
2902 efficiencysigTOFPIDD[0]->SetTitle("");
2903 efficiencysigTOFPIDD[0]->SetStats(0);
2904 efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2905 efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2906 efficiencysigTOFPIDD[0]->SetLineColor(2);
2907 efficiencysigTOFPIDD[0]->Draw();
2910 efficiencyTOFPIDD[0]->SetTitle("");
2911 efficiencyTOFPIDD[0]->SetStats(0);
2912 efficiencyTOFPIDD[0]->SetMarkerStyle(24);
2913 efficiencyTOFPIDD[0]->SetMarkerColor(4);
2914 efficiencyTOFPIDD[0]->SetLineColor(4);
2915 efficiencyTOFPIDD[0]->Draw("same");
2918 efficiencysigesdTOFPIDD[0]->SetTitle("");
2919 efficiencysigesdTOFPIDD[0]->SetStats(0);
2920 efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2921 efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2922 efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2923 efficiencysigesdTOFPIDD[0]->Draw("same");
2926 efficiencyesdTOFPIDD[0]->SetTitle("");
2927 efficiencyesdTOFPIDD[0]->SetStats(0);
2928 efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2929 efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2930 efficiencyesdTOFPIDD[0]->SetLineColor(1);
2931 efficiencyesdTOFPIDD[0]->Draw("same");
2934 if(fEfficiencyTOFPIDD[0]){
2935 fEfficiencyTOFPIDD[0]->SetLineColor(2);
2936 fEfficiencyTOFPIDD[0]->Draw("same");
2939 if(fEfficiencyesdTOFPIDD[0]){
2940 fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2941 fEfficiencyesdTOFPIDD[0]->Draw("same");
2944 TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2945 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2946 legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2947 legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2948 legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2949 legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2950 legtofeff->Draw("same");
2953 TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
2954 cefficiencyParamIP->cd();
2955 gStyle->SetOptStat(0);
2957 // IP cut efficiencies
2958 for(int icentr=0; icentr<loopcentr; icentr++) {
2960 AliCFContainer *charmCombined = 0x0;
2961 AliCFContainer *beautyCombined = 0x0;
2962 AliCFContainer *beautyCombinedesd = 0x0;
2964 printf("centrality printing %i \n",icentr);
2966 source = 0; //charm enhenced
2967 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2968 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2969 AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2970 efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2972 charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
2974 source = 1; //beauty enhenced
2975 if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2976 else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2977 AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2978 efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2980 beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
2982 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2983 else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2984 AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
2985 efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2987 beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
2989 source = 0; //charm mb
2990 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2991 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2992 AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2993 efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2995 charmCombined->Add(mccontainermcD);
2996 AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
2997 efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2999 source = 1; //beauty mb
3000 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3001 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3002 AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
3003 efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3005 beautyCombined->Add(mccontainermcD);
3006 AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
3007 efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3009 if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
3010 else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3011 AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
3012 efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3014 beautyCombinedesd->Add(mccontaineresdD);
3015 AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
3016 efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3018 source = 2; //conversion mb
3019 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3020 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3021 AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
3022 efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3024 source = 3; //non HFE except for the conversion mb
3025 if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3026 else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3027 AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
3028 efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3030 if(fIPEffCombinedSamples){
3031 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
3032 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
3033 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
3036 fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
3037 fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
3038 fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
3040 fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
3041 fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
3042 fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
3043 fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
3048 AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
3049 fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
3051 AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
3052 fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
3055 for(int icentr=0; icentr<loopcentr; icentr++) {
3056 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3057 fipfit->SetLineColor(2);
3058 fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
3059 fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3060 if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
3061 else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
3063 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3064 fipfit->SetLineColor(6);
3065 fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
3066 fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
3067 if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
3068 else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
3070 fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3071 fipfit->SetLineColor(1);
3072 fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
3073 fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3074 fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
3076 if(fIPParameterizedEff){
3077 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3078 fipfitnonhfe->SetLineColor(3);
3079 fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
3080 fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3081 fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
3083 fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3084 fipfitnonhfe->SetLineColor(4);
3085 fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
3086 fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3087 fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
3091 // draw (for PbPb, only 1st bin)
3092 fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
3093 fEfficiencyCharmSigD[0]->SetMarkerColor(1);
3094 fEfficiencyCharmSigD[0]->SetLineColor(1);
3095 fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
3096 fEfficiencyBeautySigD[0]->SetMarkerColor(2);
3097 fEfficiencyBeautySigD[0]->SetLineColor(2);
3098 fEfficiencyBeautySigesdD[0]->SetStats(0);
3099 fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
3100 fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
3101 fEfficiencyBeautySigesdD[0]->SetLineColor(6);
3102 fCharmEff[0]->SetMarkerStyle(24);
3103 fCharmEff[0]->SetMarkerColor(1);
3104 fCharmEff[0]->SetLineColor(1);
3105 fBeautyEff[0]->SetMarkerStyle(24);
3106 fBeautyEff[0]->SetMarkerColor(2);
3107 fBeautyEff[0]->SetLineColor(2);
3108 fConversionEff[0]->SetMarkerStyle(24);
3109 fConversionEff[0]->SetMarkerColor(3);
3110 fConversionEff[0]->SetLineColor(3);
3111 fNonHFEEff[0]->SetMarkerStyle(24);
3112 fNonHFEEff[0]->SetMarkerColor(4);
3113 fNonHFEEff[0]->SetLineColor(4);
3115 fEfficiencyCharmSigD[0]->Draw();
3116 fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
3117 fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
3119 fEfficiencyBeautySigD[0]->Draw("same");
3120 fEfficiencyBeautySigesdD[0]->Draw("same");
3121 //fCharmEff[0]->Draw("same");
3122 //fBeautyEff[0]->Draw("same");
3125 fConversionEffbgc->SetMarkerStyle(25);
3126 fConversionEffbgc->SetMarkerColor(3);
3127 fConversionEffbgc->SetLineColor(3);
3128 fNonHFEEffbgc->SetMarkerStyle(25);
3129 fNonHFEEffbgc->SetMarkerColor(4);
3130 fNonHFEEffbgc->SetLineColor(4);
3131 fConversionEffbgc->Draw("same");
3132 fNonHFEEffbgc->Draw("same");
3135 fConversionEff[0]->Draw("same");
3136 fNonHFEEff[0]->Draw("same");
3138 if(fEfficiencyIPBeautyD[0])
3139 fEfficiencyIPBeautyD[0]->Draw("same");
3140 if(fEfficiencyIPBeautyesdD[0])
3141 fEfficiencyIPBeautyesdD[0]->Draw("same");
3142 if( fEfficiencyIPCharmD[0])
3143 fEfficiencyIPCharmD[0]->Draw("same");
3144 if(fIPParameterizedEff){
3145 fEfficiencyIPConversionD[0]->Draw("same");
3146 fEfficiencyIPNonhfeD[0]->Draw("same");
3148 TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
3149 legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
3150 legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
3151 legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
3152 legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
3153 legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
3154 legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
3155 //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
3156 //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
3157 legipeff->Draw("same");
3159 //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
3162 //____________________________________________________________________________
3163 THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
3165 // Return beauty electron IP cut efficiency
3168 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3169 const Int_t nCentralitybinning=11;//number of centrality bins
3170 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
3171 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3173 Int_t nDim=1; //dimensions of the efficiency weighting grid
3177 nDim=2; //dimensions of the efficiency weighting grid
3179 Int_t nBin[1] = {nPtbinning1};
3180 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3184 if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
3185 else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3187 for(Int_t idim = 0; idim < nDim; idim++)
3189 if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
3192 ipcut->SetBinEdges(0, kCentralityRange);
3193 ipcut->SetBinEdges(1, kPtRange);
3197 Double_t weight[nCentralitybinning];
3198 Double_t weightErr[nCentralitybinning];
3199 Double_t contents[2];
3201 for(Int_t a=0;a<11;a++)
3208 Int_t looppt=nBin[0];
3213 loopcentr=nBinPbPb[0];
3217 for(int icentr=0; icentr<loopcentr; icentr++)
3219 for(int i=0; i<looppt; i++)
3221 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3223 if(fEfficiencyIPBeautyD[icentr]){
3224 weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
3225 weightErr[icentr] = 0;
3228 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3229 weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
3230 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3234 if(fEfficiencyIPBeautyesdD[icentr]){
3235 weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
3236 weightErr[icentr] = 0;
3239 printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3240 weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
3241 weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3255 ipcut->Fill(contents,weight[icentr]);
3256 ipcut->SetBinError(ibin,weightErr[icentr]);
3260 Int_t nDimSparse = ipcut->GetNdimensions();
3261 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3262 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3264 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3265 binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
3266 nBins *= binsvar[iVar];
3269 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3270 // loop that sets 0 error in each bin
3271 for (Long_t iBin=0; iBin<nBins; iBin++) {
3272 Long_t bintmp = iBin ;
3273 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3274 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3275 bintmp /= binsvar[iVar] ;
3277 //ipcut->SetBinError(binfill,0.); // put 0 everywhere
3286 //____________________________________________________________________________
3287 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3289 // Return PID x IP cut efficiency
3291 const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3292 const Int_t nCentralitybinning=11;//number of centrality bins
3293 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
3294 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3296 Int_t nDim=1; //dimensions of the efficiency weighting grid
3300 nDim=2; //dimensions of the efficiency weighting grid
3302 Int_t nBin[1] = {nPtbinning1};
3303 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3306 if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3307 else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3308 for(Int_t idim = 0; idim < nDim; idim++)
3311 if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3314 pideff->SetBinEdges(0, kCentralityRange);
3315 pideff->SetBinEdges(1, kPtRange);
3320 Double_t weight[nCentralitybinning];
3321 Double_t weightErr[nCentralitybinning];
3322 Double_t contents[2];
3324 for(Int_t a=0;a<nCentralitybinning;a++)
3330 Int_t looppt=nBin[0];
3335 loopcentr=nBinPbPb[0];
3338 for(int icentr=0; icentr<loopcentr; icentr++)
3340 Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3341 for(int i=0; i<looppt; i++)
3343 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3345 Double_t tofpideff = 0.;
3346 Double_t tofpideffesd = 0.;
3347 if(fEfficiencyTOFPIDD[icentr])
3348 tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
3350 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3352 if(fEfficiencyesdTOFPIDD[icentr])
3353 tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3355 printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3358 //tof pid eff x tpc pid eff x ip cut eff
3359 if(fIPParameterizedEff){
3361 if(fEfficiencyIPCharmD[icentr]){
3362 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3363 weightErr[icentr] = 0;
3366 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3367 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3368 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3371 else if(source==2) {
3372 if(fEfficiencyIPConversionD[icentr]){
3373 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
3374 weightErr[icentr] = 0;
3377 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
3378 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3379 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3382 else if(source==3) {
3383 if(fEfficiencyIPNonhfeD[icentr]){
3384 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
3385 weightErr[icentr] = 0;
3388 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
3389 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3390 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3396 if(fEfficiencyIPCharmD[icentr]){
3397 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3398 weightErr[icentr] = 0;
3401 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3402 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3403 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3408 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3409 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3412 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3413 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3418 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3419 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3422 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3423 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3439 pideff->Fill(contents,weight[icentr]);
3440 pideff->SetBinError(ibin,weightErr[icentr]);
3444 Int_t nDimSparse = pideff->GetNdimensions();
3445 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3446 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3448 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3449 binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3450 nBins *= binsvar[iVar];
3453 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3454 // loop that sets 0 error in each bin
3455 for (Long_t iBin=0; iBin<nBins; iBin++) {
3456 Long_t bintmp = iBin ;
3457 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3458 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3459 bintmp /= binsvar[iVar] ;
3470 //__________________________________________________________________________
3471 AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
3473 // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
3487 const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
3488 const Int_t nCentralitybinning=11;//number of centrality bins
3489 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};
3491 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3492 Int_t nBin[1] = {nPtbinning1};
3493 Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3495 AliCFDataGrid *rawBeautyContainer;
3496 if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
3497 else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
3498 // printf("number of bins= %d\n",bins[0]);
3503 THnSparseF *brawspectra;
3504 if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
3505 else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
3506 if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
3509 // brawspectra->SetBinEdges(0, centralityBins);
3510 brawspectra->SetBinEdges(0, kCentralityRange);
3511 brawspectra->SetBinEdges(1, kPtRange);
3515 Double_t yields= 0.;
3516 Double_t valuesb[2];
3518 //Int_t looppt=nBin[0];
3522 loopcentr=nBinPbPb[0];
3525 for(int icentr=0; icentr<loopcentr; icentr++)
3528 for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
3529 pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3531 yields = fBSpectrum2ndMethod->GetBinContent(i+1);
3538 if(fBeamType==0) valuesb[0]=pt[0];
3539 brawspectra->Fill(valuesb,yields);
3545 Int_t nDimSparse = brawspectra->GetNdimensions();
3546 Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3547 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3549 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3550 binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
3551 nBins *= binsvar[iVar];
3554 Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3555 // loop that sets 0 error in each bin
3556 for (Long_t iBin=0; iBin<nBins; iBin++) {
3557 Long_t bintmp = iBin ;
3558 for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3559 binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3560 bintmp /= binsvar[iVar] ;
3562 brawspectra->SetBinError(binfill,0.); // put 0 everywhere
3566 rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
3567 TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
3570 fBSpectrum2ndMethod->SetMarkerStyle(24);
3571 fBSpectrum2ndMethod->Draw("p");
3572 hRawBeautySpectra->SetMarkerStyle(25);
3573 hRawBeautySpectra->Draw("samep");
3578 return rawBeautyContainer;
3581 //__________________________________________________________________________
3582 void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
3584 // Calculate non HFE sys
3591 Double_t evtnorm[1] = {0.0};
3592 if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
3594 AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
3595 AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
3597 AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
3598 AliCFDataGrid *bgNonHFEGrid[kBgLevels];
3599 AliCFDataGrid *bgConvGrid[kBgLevels];
3601 Int_t stepbackground = 3;
3602 Int_t* bins=new Int_t[1];
3603 const Char_t *bgBase[2] = {"pi0","eta"};
3605 bins[0]=fConversionEff[centrality]->GetNbinsX();
3607 AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
3608 AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
3610 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
3611 for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
3612 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);
3613 weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
3614 convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
3616 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);
3617 weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
3618 nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
3621 bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
3622 for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
3623 bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
3626 bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
3628 bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
3629 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)
3630 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
3633 bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
3635 bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
3636 bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
3638 bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
3639 bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
3644 //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)
3645 AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
3646 TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
3647 for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
3648 bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
3649 bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
3650 bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
3651 bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
3653 //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
3655 hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
3656 hBaseErrors[iErr][0]->Scale(-1.);
3657 hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
3658 hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
3663 //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
3664 TH1D *hSpeciesErrors[kElecBgSources-1];
3665 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3666 if(fEtaSyst && (iSource == 1))continue;
3667 hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
3668 TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
3669 hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
3670 hSpeciesErrors[iSource-1]->Scale(0.3);
3673 //Int_t firstBgSource = 0;//if eta systematics are not from scaling
3674 //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
3675 TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
3676 TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
3677 TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
3679 TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
3680 TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
3682 for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
3683 Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
3684 pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
3685 pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
3687 etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
3688 etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
3690 else{ etaErrLow = etaErrUp = 0.;}
3692 Double_t sqrsumErrs= 0;
3693 for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3694 if(fEtaSyst && (iSource == 1))continue;
3695 Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
3696 sqrsumErrs+=(scalingErr*scalingErr);
3698 for(Int_t iErr = 0; iErr < 2; iErr++){
3699 for(Int_t iLevel = 0; iLevel < 2; iLevel++){
3700 hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
3704 hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
3705 hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
3706 hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
3708 hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
3709 hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
3713 TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
3715 cPiErrors->SetLogx();
3716 cPiErrors->SetLogy();
3717 hBaseErrors[0][0]->Draw();
3718 //hBaseErrors[0][1]->SetMarkerColor(kBlack);
3719 //hBaseErrors[0][1]->SetLineColor(kBlack);
3720 //hBaseErrors[0][1]->Draw("SAME");
3722 hBaseErrors[1][0]->Draw("SAME");
3723 hBaseErrors[1][0]->SetMarkerColor(kBlack);
3724 hBaseErrors[1][0]->SetLineColor(kBlack);
3725 //hBaseErrors[1][1]->SetMarkerColor(13);
3726 //hBaseErrors[1][1]->SetLineColor(13);
3727 //hBaseErrors[1][1]->Draw("SAME");
3729 //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
3730 //hOverallBinScaledErrsUp->SetLineColor(kBlue);
3731 //hOverallBinScaledErrsUp->Draw("SAME");
3732 hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
3733 hOverallBinScaledErrsLow->SetLineColor(kGreen);
3734 hOverallBinScaledErrsLow->Draw("SAME");
3735 hScalingErrors->SetLineColor(kBlue);
3736 hScalingErrors->Draw("SAME");
3738 TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
3739 lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
3740 //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
3742 lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
3743 //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
3745 lPiErr->AddEntry(hScalingErrors, "scaling error");
3746 lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
3747 //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
3748 lPiErr->Draw("SAME");
3751 TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
3752 TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
3753 hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
3754 hLowSystScaled->Scale(evtnorm[0]);
3755 TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
3756 TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
3757 //histograms to be normalized to TGraphErrors
3758 CorrectFromTheWidth(hNormAllSystErrUp);
3759 CorrectFromTheWidth(hNormAllSystErrLow);
3761 TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
3763 cNormOvErrs->SetLogx();
3764 cNormOvErrs->SetLogy();
3766 TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
3767 TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
3768 gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
3769 gOverallSystErrUp->SetMarkerColor(kBlack);
3770 gOverallSystErrUp->SetLineColor(kBlack);
3771 gOverallSystErrLow->SetMarkerColor(kRed);
3772 gOverallSystErrLow->SetLineColor(kRed);
3773 gOverallSystErrUp->Draw("AP");
3774 gOverallSystErrLow->Draw("PSAME");
3775 TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
3776 lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
3777 lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
3778 lAllSys->Draw("same");
3781 AliCFDataGrid *bgYieldGrid;
3783 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.
3785 bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
3787 TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
3788 TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
3789 hRelErrUp->Divide(hBgYield);
3790 TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
3791 hRelErrLow->Divide(hBgYield);
3793 TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
3795 cRelErrs->SetLogx();
3796 hRelErrUp->SetTitle("Relative error of non-HFE background yield");
3798 hRelErrLow->SetLineColor(kBlack);
3799 hRelErrLow->Draw("SAME");
3801 TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
3802 lRel->AddEntry(hRelErrUp, "upper");
3803 lRel->AddEntry(hRelErrLow, "lower");
3806 //CorrectFromTheWidth(hBgYield);
3807 //hBgYield->Scale(evtnorm[0]);
3810 //write histograms/TGraphs to file
3811 TFile *output = new TFile("systHists.root","recreate");
3813 hBgYield->SetName("hBgYield");
3815 hRelErrUp->SetName("hRelErrUp");
3817 hRelErrLow->SetName("hRelErrLow");
3818 hRelErrLow->Write();
3819 hUpSystScaled->SetName("hOverallSystErrUp");
3820 hUpSystScaled->Write();
3821 hLowSystScaled->SetName("hOverallSystErrLow");
3822 hLowSystScaled->Write();
3823 gOverallSystErrUp->SetName("gOverallSystErrUp");
3824 gOverallSystErrUp->Write();
3825 gOverallSystErrLow->SetName("gOverallSystErrLow");
3826 gOverallSystErrLow->Write();
3833 //____________________________________________________________
3834 void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
3837 // Unfold backgrounds to check its sanity
3840 AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
3841 //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
3843 AliError("MC Container not available");
3848 AliError("No Correlation map available");
3853 AliCFDataGrid *dataGrid = 0x0;
3854 dataGrid = bgsubpectrum;
3857 AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
3858 THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
3861 AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
3862 efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
3864 // Unfold background spectra
3866 if(fBeamType==0)nDim = 1;
3867 if(fBeamType==1)nDim = 2;
3868 AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
3869 if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
3870 unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
3871 if(fSetSmoothing) unfolding.UseSmoothing();
3875 THnSparse* result = unfolding.GetUnfolded();
3876 TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
3881 result->GetAxis(0)->SetRange(1,1);
3882 TH1D* htest1=(TH1D*)result->Projection(0);
3885 result->GetAxis(0)->SetRange(1,1);
3886 TH1D* htest2=(TH1D*)result->Projection(1);
3889 result->GetAxis(0)->SetRange(6,6);
3890 TH1D* htest3=(TH1D*)result->Projection(0);
3893 result->GetAxis(0)->SetRange(6,6);
3894 TH1D* htest4=(TH1D*)result->Projection(1);
3903 TGraphErrors* unfoldedbgspectrumD = Normalize(result);
3904 if(!unfoldedbgspectrumD) {
3905 AliError("Unfolded background spectrum doesn't exist");
3908 TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
3909 if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
3914 result->GetAxis(0)->SetRange(centr,centr);
3915 unfoldedbgspectrumD = Normalize(result,centr-1);
3916 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
3918 result->GetAxis(0)->SetRange(centr,centr);
3919 unfoldedbgspectrumD = Normalize(result,centr-1);
3920 unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");