/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ // // Class for spectrum correction // Subtraction of hadronic background, Unfolding of the data and // Renormalization done here // The following containers have to be set: // - Correction framework container for real data // - Correction framework container for MC (Efficiency Map) // - Correction framework container for background coming from data // - Correction framework container for background coming from MC // // Author: // Raphaelle Bailhache // Markus Fasel // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "AliPID.h" #include "AliCFContainer.h" #include "AliCFDataGrid.h" #include "AliCFEffGrid.h" #include "AliCFGridSparse.h" #include "AliCFUnfolding.h" #include "AliLog.h" #include "AliHFEspectrum.h" #include "AliHFEcuts.h" #include "AliHFEcontainer.h" #include "AliHFEtools.h" ClassImp(AliHFEspectrum) //____________________________________________________________ AliHFEspectrum::AliHFEspectrum(const char *name): TNamed(name, ""), fCFContainers(new TObjArray(kDataContainerV0+1)), fTemporaryObjects(NULL), fCorrelation(NULL), fBackground(NULL), fEfficiencyFunction(NULL), fWeightCharm(NULL), fInclusiveSpectrum(kTRUE), fDumpToFile(kFALSE), fEtaSelected(kFALSE), fUnSetCorrelatedErrors(kTRUE), fSetSmoothing(kFALSE), fIPanaHadronBgSubtract(kFALSE), fIPanaCharmBgSubtract(kFALSE), fIPanaConversionBgSubtract(kFALSE), fIPanaNonHFEBgSubtract(kFALSE), fIPParameterizedEff(kFALSE), fNonHFEsyst(0), fBeauty2ndMethod(kFALSE), fIPEffCombinedSamples(kTRUE), fNbDimensions(1), fNMCEvents(0), fStepMC(-1), fStepTrue(-1), fStepData(-1), fStepBeforeCutsV0(-1), fStepAfterCutsV0(-1), fStepGuessedUnfolding(-1), fNumberOfIterations(5), fNRandomIter(0), fChargeChoosen(kAllCharge), fNCentralityBinAtTheEnd(0), fTestCentralityLow(-1), fTestCentralityHigh(-1), fFillMoreCorrelationMatrix(kFALSE), fHadronEffbyIPcut(NULL), fConversionEffbgc(NULL), fNonHFEEffbgc(NULL), fBSpectrum2ndMethod(NULL), fkBeauty2ndMethodfilename(""), fBeamType(0), fEtaSyst(kTRUE), fDebugLevel(0), fWriteToFile(kFALSE), fUnfoldBG(kFALSE) { // // Default constructor // for(Int_t k = 0; k < 20; k++){ fNEvents[k] = 0; fNMCbgEvents[k] = 0; fLowBoundaryCentralityBinAtTheEnd[k] = 0; fHighBoundaryCentralityBinAtTheEnd[k] = 0; if(kClear(); delete fTemporaryObjects; } } //____________________________________________________________ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){ // // Init what we need for the correction: // // Raw spectrum, hadron contamination // MC efficiency maps, correlation matrix // V0 efficiency if wanted // // This for a given dimension. // If no inclusive spectrum, then take only efficiency map for beauty electron // and the appropriate correlation matrix // if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod(); Int_t kNdim = 3; Int_t kNcentr =1; //Int_t ptpr =0; if(fBeamType==0) kNdim=3; if(fBeamType==1) { kNdim=4; kNcentr=11; //ptpr=1; } Int_t dims[kNdim]; // Get the requested format if(fBeamType==0){ // pp analysis switch(fNbDimensions){ case 1: dims[0] = 0; break; case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i; break; case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i; break; default: AliError("Container with this number of dimensions not foreseen (yet)"); return kFALSE; }; } if(fBeamType==1){ // PbPb analysis; centrality as first dimension Int_t nbDimensions = fNbDimensions; fNbDimensions = fNbDimensions + 1; switch(nbDimensions){ case 1: dims[0] = 5; dims[1] = 0; break; case 2: dims[0] = 5; for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i; break; case 3: dims[0] = 5; for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i; break; default: AliError("Container with this number of dimensions not foreseen (yet)"); return kFALSE; }; } // Data container: raw spectrum + hadron contamination AliCFContainer *datacontainer = 0x0; if(fInclusiveSpectrum) { datacontainer = datahfecontainer->GetCFContainer("recTrackContReco"); } else{ datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco"); } AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground"); if((!datacontainer) || (!contaminationcontainer)) return kFALSE; AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE; SetContainer(datacontainerD,AliHFEspectrum::kDataContainer); SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData); // MC container: ESD/MC efficiency maps + MC/MC efficiency maps AliCFContainer *mccontaineresd = 0x0; AliCFContainer *mccontaineresdbg = 0x0; AliCFContainer *mccontainermc = 0x0; AliCFContainer *mccontainermcbg = 0x0; AliCFContainer *nonHFEweightedContainer = 0x0; AliCFContainer *convweightedContainer = 0x0; AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers if(fInclusiveSpectrum) { mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco"); mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC"); } else { mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco"); mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco"); mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC"); mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC"); if(fNonHFEsyst){ const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"}; const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"}; for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){ for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel])); convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel])); for(Int_t icentr=0;icentrGetCFContainer("mesonElecs"); convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs"); if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE; //} } if((!mccontaineresd) || (!mccontainermc)) return kFALSE; Int_t source = -1; if(!fInclusiveSpectrum) source = 1; //beauty AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE; SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC); SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD); // set charm, nonHFE container to estimate BG if(!fInclusiveSpectrum){ source = 0; //charm mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen); SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC); //if(!fNonHFEsyst){ AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen); SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD); AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen); SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD); //} SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims); } // MC container: correlation matrix THnSparseF *mccorrelation = 0x0; if(fInclusiveSpectrum) { if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID"); else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID"); else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID"); else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID"); else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); } else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE"); if(!mccorrelation) return kFALSE; Int_t testCentralityLow = fTestCentralityLow; Int_t testCentralityHigh = fTestCentralityHigh; if(fFillMoreCorrelationMatrix) { testCentralityLow = fTestCentralityLow-1; testCentralityHigh = fTestCentralityHigh+1; } THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh); if(!mccorrelationD) { printf("No correlation\n"); return kFALSE; } SetCorrelation(mccorrelationD); // V0 container Electron, pt eta phi if(v0hfecontainer) { AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco"); if(!containerV0) return kFALSE; AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); if(!containerV0Electron) return kFALSE; SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0); } if(fDebugLevel>0){ AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone(); contaminationspectrum->SetName("contaminationspectrum"); TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700); ccontaminationspectrum->Divide(3,1); ccontaminationspectrum->cd(1); TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0); TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0); TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2); contaminationspectrum2dpteta->SetStats(0); contaminationspectrum2dpteta->SetTitle(""); contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta"); contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]"); contaminationspectrum2dptphi->SetStats(0); contaminationspectrum2dptphi->SetTitle(""); contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]"); contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]"); contaminationspectrum2detaphi->SetStats(0); contaminationspectrum2detaphi->SetTitle(""); contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta"); contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]"); contaminationspectrum2dptphi->Draw("colz"); ccontaminationspectrum->cd(2); contaminationspectrum2dpteta->Draw("colz"); ccontaminationspectrum->cd(3); contaminationspectrum2detaphi->Draw("colz"); TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700); ccorrelation->cd(1); if(mccorrelationD) { if(fBeamType==0){ TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0); ptcorrelation->Draw("colz"); } if(fBeamType==1){ TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1); ptcorrelation->Draw("colz"); } } if(fWriteToFile){ ccontaminationspectrum->SaveAs("contaminationspectrum.eps"); ccorrelation->SaveAs("correlationMatrix.eps"); } } TFile *file = TFile::Open("tests.root","recreate"); datacontainerD->Write("data"); mccontainermcD->Write("mcefficiency"); mccorrelationD->Write("correlationmatrix"); file->Close(); return kTRUE; } //____________________________________________________________ void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){ // // get spectrum for beauty 2nd method // // TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename); fBSpectrum2ndMethod = new TH1D(*static_cast(inputfile2ndmethod->Get("BSpectrum"))); } //____________________________________________________________ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){ // // Correct the spectrum for efficiency and unfolding // with both method and compare // gStyle->SetPalette(1); gStyle->SetOptStat(1111); gStyle->SetPadBorderMode(0); gStyle->SetCanvasColor(10); gStyle->SetPadLeftMargin(0.13); gStyle->SetPadRightMargin(0.13); printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0); /////////////////////////// // Check initialization /////////////////////////// if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){ AliInfo("You have to init before"); return kFALSE; } if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) { AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect"); return kFALSE; } SetNumberOfIteration(10); SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack); AliCFDataGrid *dataGridAfterFirstSteps = 0x0; ////////////////////////////////// // Subtract hadron background ///////////////////////////////// AliCFDataGrid *dataspectrumaftersubstraction = 0x0; if(subtractcontamination) { dataspectrumaftersubstraction = SubtractBackground(kTRUE); dataGridAfterFirstSteps = dataspectrumaftersubstraction; } printf("cloning spectrum\n"); AliCFDataGrid *rawsave(NULL); if(dataspectrumaftersubstraction){ rawsave = (AliCFDataGrid *)dataspectrumaftersubstraction->Clone("rawdata"); } else { AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); } rawsave = new AliCFDataGrid("rawsave", "raw spectrum after subtraction",*dataContainer, fStepData); } printf("cloned: %p\n", rawsave); //////////////////////////////////////////////// // Correct for TPC efficiency from V0 /////////////////////////////////////////////// AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0; AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0); if(dataContainerV0){printf("Got the V0 container\n"); dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction); dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection; } ////////////////////////////////////////////////////////////////////////////// // Correct for efficiency parametrized (if TPC efficiency is parametrized) ///////////////////////////////////////////////////////////////////////////// AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0; if(fEfficiencyFunction){ dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps); dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection; } /////////////// // Unfold ////////////// TList *listunfolded = Unfold(dataGridAfterFirstSteps); if(!listunfolded){ printf("Unfolded failed\n"); return kFALSE; } THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0); THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1); if(!correctedspectrum){ AliError("No corrected spectrum\n"); return kFALSE; } if(!residualspectrum){ AliError("No residul spectrum\n"); return kFALSE; } ///////////////////// // Simply correct //////////////////// AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps); ////////// // Plot ////////// if(fDebugLevel > 0.0) { Int_t ptpr = 0; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700); ccorrected->Divide(2,1); ccorrected->cd(1); gPad->SetLogy(); TGraphErrors* correctedspectrumD = Normalize(correctedspectrum); correctedspectrumD->SetTitle(""); correctedspectrumD->GetYaxis()->SetTitleOffset(1.5); correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); correctedspectrumD->SetMarkerStyle(26); correctedspectrumD->SetMarkerColor(kBlue); correctedspectrumD->SetLineColor(kBlue); correctedspectrumD->Draw("AP"); TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection); alltogetherspectrumD->SetTitle(""); alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5); alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); alltogetherspectrumD->SetMarkerStyle(25); alltogetherspectrumD->SetMarkerColor(kBlack); alltogetherspectrumD->SetLineColor(kBlack); alltogetherspectrumD->Draw("P"); TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89); legcorrected->AddEntry(correctedspectrumD,"Corrected","p"); legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p"); legcorrected->Draw("same"); ccorrected->cd(2); TH1D *correctedTH1D = correctedspectrum->Projection(ptpr); TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr); TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone(); ratiocorrected->SetName("ratiocorrected"); ratiocorrected->SetTitle(""); ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected"); ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1); ratiocorrected->SetStats(0); ratiocorrected->Draw(); if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps"); //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd]; //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd]; //TH1D correctedspectrac[fNCentralityBinAtTheEnd]; //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd]; TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd]; TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd]; TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd]; TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd]; if(fBeamType==1) { TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700); ccorrectedallspectra->Divide(2,1); TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89); TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89); THnSparseF* sparsesu = (THnSparseF *) correctedspectrum; TAxis *cenaxisa = sparsesu->GetAxis(0); THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid(); TAxis *cenaxisb = sparsed->GetAxis(0); Int_t nbbin = cenaxisb->GetNbins(); Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29}; Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20}; for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){ TString titlee("corrected_centrality_bin_"); titlee += "["; titlee += fLowBoundaryCentralityBinAtTheEnd[binc]; titlee += "_"; titlee += fHighBoundaryCentralityBinAtTheEnd[binc]; titlee += "["; TString titleec("corrected_check_projection_bin_"); titleec += "["; titleec += fLowBoundaryCentralityBinAtTheEnd[binc]; titleec += "_"; titleec += fHighBoundaryCentralityBinAtTheEnd[binc]; titleec += "["; TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_"); titleenameunotnormalized += "["; titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc]; titleenameunotnormalized += "_"; titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc]; titleenameunotnormalized += "["; TString titleenameunormalized("Unfolded_normalized_centrality_bin_"); titleenameunormalized += "["; titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc]; titleenameunormalized += "_"; titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc]; titleenameunormalized += "["; TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_"); titleenamednotnormalized += "["; titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc]; titleenamednotnormalized += "_"; titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc]; titleenamednotnormalized += "["; TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_"); titleenamednormalized += "["; titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc]; titleenamednormalized += "_"; titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc]; titleenamednormalized += "["; Int_t nbEvents = 0; for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) { printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k); nbEvents += fNEvents[k]; } Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1); Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]); printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega); Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1); Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]); printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb); cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]); cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]); TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700); ccorrectedcheck->cd(1); TH1D *aftersuc = (TH1D *) sparsesu->Projection(0); TH1D *aftersdc = (TH1D *) sparsed->Projection(0); aftersuc->Draw(); aftersdc->Draw("same"); TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700); ccorrectede->Divide(2,1); ccorrectede->cd(1); gPad->SetLogy(); TH1D *aftersu = (TH1D *) sparsesu->Projection(1); CorrectFromTheWidth(aftersu); aftersu->SetName((const char*)titleenameunotnormalized); unfoldingspectrac[binc] = *aftersu; ccorrectede->cd(1); TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents); aftersun->SetTitle(""); aftersun->GetYaxis()->SetTitleOffset(1.5); aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0); aftersun->SetMarkerStyle(26); aftersun->SetMarkerColor(kBlue); aftersun->SetLineColor(kBlue); aftersun->Draw("AP"); aftersun->SetName((const char*)titleenameunormalized); unfoldingspectracn[binc] = *aftersun; ccorrectede->cd(1); TH1D *aftersd = (TH1D *) sparsed->Projection(1); CorrectFromTheWidth(aftersd); aftersd->SetName((const char*)titleenamednotnormalized); correctedspectrac[binc] = *aftersd; ccorrectede->cd(1); TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents); aftersdn->SetTitle(""); aftersdn->GetYaxis()->SetTitleOffset(1.5); aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0); aftersdn->SetMarkerStyle(25); aftersdn->SetMarkerColor(kBlack); aftersdn->SetLineColor(kBlack); aftersdn->Draw("P"); aftersdn->SetName((const char*)titleenamednormalized); correctedspectracn[binc] = *aftersdn; TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89); legcorrectedud->AddEntry(aftersun,"Corrected","p"); legcorrectedud->AddEntry(aftersdn,"Alltogether","p"); legcorrectedud->Draw("same"); ccorrectedallspectra->cd(1); gPad->SetLogy(); TH1D *aftersunn = (TH1D *) aftersun->Clone(); aftersunn->SetMarkerStyle(stylee[binc]); aftersunn->SetMarkerColor(colorr[binc]); if(binc==0) aftersunn->Draw("AP"); else aftersunn->Draw("P"); legtotal->AddEntry(aftersunn,(const char*) titlee,"p"); ccorrectedallspectra->cd(2); gPad->SetLogy(); TH1D *aftersdnn = (TH1D *) aftersdn->Clone(); aftersdnn->SetMarkerStyle(stylee[binc]); aftersdnn->SetMarkerColor(colorr[binc]); if(binc==0) aftersdnn->Draw("AP"); else aftersdnn->Draw("P"); legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p"); ccorrectede->cd(2); TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone(); TString titleee("ratiocorrected_bin_"); titleee += binc; ratiocorrectedbinc->SetName((const char*) titleee); ratiocorrectedbinc->SetTitle(""); ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected"); ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratiocorrectedbinc->Divide(aftersu,aftersd,1,1); ratiocorrectedbinc->SetStats(0); ratiocorrectedbinc->Draw(); } ccorrectedallspectra->cd(1); legtotal->Draw("same"); ccorrectedallspectra->cd(2); legtotalg->Draw("same"); cenaxisa->SetRange(0,nbbin); cenaxisb->SetRange(0,nbbin); if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps"); } // Dump to file if needed if(fDumpToFile) { TFile *out = new TFile("finalSpectrum.root","recreate"); correctedspectrumD->SetName("UnfoldingCorrectedSpectrum"); correctedspectrumD->Write(); alltogetherspectrumD->SetName("AlltogetherSpectrum"); alltogetherspectrumD->Write(); ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum"); ratiocorrected->Write(); correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum"); correctedspectrum->Write(); alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum"); alltogetherCorrection->Write(); rawsave->Write(); for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){ unfoldingspectrac[binc].Write(); unfoldingspectracn[binc].Write(); correctedspectrac[binc].Write(); correctedspectracn[binc].Write(); } out->Close(); delete out; } if (unfoldingspectrac) delete[] unfoldingspectrac; if (unfoldingspectracn) delete[] unfoldingspectracn; if (correctedspectrac) delete[] correctedspectrac; if (correctedspectracn) delete[] correctedspectracn; } return kTRUE; } //____________________________________________________________ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ // // Correct the spectrum for efficiency and unfolding for beauty analysis // with both method and compare // gStyle->SetPalette(1); gStyle->SetOptStat(1111); gStyle->SetPadBorderMode(0); gStyle->SetCanvasColor(10); gStyle->SetPadLeftMargin(0.13); gStyle->SetPadRightMargin(0.13); /////////////////////////// // Check initialization /////////////////////////// if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){ AliInfo("You have to init before"); return kFALSE; } if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) { AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect"); return kFALSE; } SetNumberOfIteration(10); SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack); AliCFDataGrid *dataGridAfterFirstSteps = 0x0; ////////////////////////////////// // Subtract hadron background ///////////////////////////////// AliCFDataGrid *dataspectrumaftersubstraction = 0x0; AliCFDataGrid *unnormalizedRawSpectrum = 0x0; TGraphErrors *gNormalizedRawSpectrum = 0x0; if(subtractcontamination) { if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE); else dataspectrumaftersubstraction = GetRawBspectra2ndMethod(); unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone(); dataGridAfterFirstSteps = dataspectrumaftersubstraction; gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum); } printf("after normalize getting IP \n"); ///////////////////////////////////////////////////////////////////////////////////////// // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds ///////////////////////////////////////////////////////////////////////////////////////// AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0; AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0; AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0); if(fEfficiencyFunction){ dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps); dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection; } else if(dataContainerV0){ dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction); dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection; } /////////////// // Unfold ////////////// TList *listunfolded = Unfold(dataGridAfterFirstSteps); if(!listunfolded){ printf("Unfolded failed\n"); return kFALSE; } THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0); THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1); if(!correctedspectrum){ AliError("No corrected spectrum\n"); return kFALSE; } if(!residualspectrum){ AliError("No residual spectrum\n"); return kFALSE; } ///////////////////// // Simply correct //////////////////// AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps); ////////// // Plot ////////// if(fDebugLevel > 0.0) { Int_t ptpr = 0; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700); ccorrected->Divide(2,1); ccorrected->cd(1); gPad->SetLogy(); TGraphErrors* correctedspectrumD = Normalize(correctedspectrum); correctedspectrumD->SetTitle(""); correctedspectrumD->GetYaxis()->SetTitleOffset(1.5); correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); correctedspectrumD->SetMarkerStyle(26); correctedspectrumD->SetMarkerColor(kBlue); correctedspectrumD->SetLineColor(kBlue); correctedspectrumD->Draw("AP"); TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection); alltogetherspectrumD->SetTitle(""); alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5); alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); alltogetherspectrumD->SetMarkerStyle(25); alltogetherspectrumD->SetMarkerColor(kBlack); alltogetherspectrumD->SetLineColor(kBlack); alltogetherspectrumD->Draw("P"); TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89); legcorrected->AddEntry(correctedspectrumD,"Corrected","p"); legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p"); legcorrected->Draw("same"); ccorrected->cd(2); TH1D *correctedTH1D = correctedspectrum->Projection(ptpr); TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr); TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone(); ratiocorrected->SetName("ratiocorrected"); ratiocorrected->SetTitle(""); ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected"); ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1); ratiocorrected->SetStats(0); ratiocorrected->Draw(); if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps"); if(fBeamType == 0){ if(fNonHFEsyst){ CalculateNonHFEsyst(0); } } // Dump to file if needed if(fDumpToFile) { // to do centrality dependent TFile *out; out = new TFile("finalSpectrum.root","recreate"); out->cd(); // correctedspectrumD->SetName("UnfoldingCorrectedSpectrum"); correctedspectrumD->Write(); alltogetherspectrumD->SetName("AlltogetherSpectrum"); alltogetherspectrumD->Write(); ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum"); ratiocorrected->Write(); // correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum"); correctedspectrum->Write(); alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum"); alltogetherCorrection->Write(); // if(unnormalizedRawSpectrum) { unnormalizedRawSpectrum->SetName("beautyAfterIP"); unnormalizedRawSpectrum->Write(); } if(gNormalizedRawSpectrum){ gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP"); gNormalizedRawSpectrum->Write(); } if(fBeamType==0) { Int_t countpp=0; fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp)); fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp)); fEfficiencyCharmSigD[countpp]->Write(); fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp)); fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp)); fEfficiencyBeautySigD[countpp]->Write(); fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp)); fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp)); fCharmEff[countpp]->Write(); fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp)); fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp)); fBeautyEff[countpp]->Write(); fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp)); fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp)); fConversionEff[countpp]->Write(); fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp)); fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp)); fNonHFEEff[countpp]->Write(); } if(fBeamType==1) { TGraphErrors* correctedspectrumDc[kCentrality]; TGraphErrors* alltogetherspectrumDc[kCentrality]; for(Int_t i=0;iGetAxis(0)->SetRange(i+1,i+1); correctedspectrumDc[i] = Normalize(correctedspectrum,i); if(correctedspectrumDc[i]){ correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i)); correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i)); correctedspectrumDc[i]->Write(); } alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1); alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i); if(alltogetherspectrumDc[i]){ alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i)); alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i)); alltogetherspectrumDc[i]->Write(); } TH1D *centrcrosscheck = correctedspectrum->Projection(0); centrcrosscheck->SetTitle(Form("centrality_%i",i)); centrcrosscheck->SetName(Form("centrality_%i",i)); centrcrosscheck->Write(); TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr); TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr); TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0); centrcrosscheck2->SetTitle(Form("centrality2_%i",i)); centrcrosscheck2->SetName(Form("centrality2_%i",i)); centrcrosscheck2->Write(); TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone(); ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1); ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i)); ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i)); ratiocorrectedc->Write(); fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i)); fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i)); fEfficiencyCharmSigD[i]->Write(); fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i)); fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i)); fEfficiencyBeautySigD[i]->Write(); fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i)); fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i)); fCharmEff[i]->Write(); fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i)); fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i)); fBeautyEff[i]->Write(); fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i)); fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i)); fConversionEff[i]->Write(); fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i)); fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i)); fNonHFEEff[i]->Write(); } } out->Close(); delete out; } } return kTRUE; } //____________________________________________________________ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ // // Apply background subtraction // Int_t ptpr = 0; Int_t nbins=1; if(fBeamType==0) { ptpr=0; nbins=1; } if(fBeamType==1) { ptpr=1; nbins=2; } // Raw spectrum AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } printf("Step data: %d\n",fStepData); AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData); AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone(); dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); // Background Estimate AliCFContainer *backgroundContainer = GetContainer(kBackgroundData); if(!backgroundContainer){ AliError("MC background container not found"); return NULL; } Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method) AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground); if(!fInclusiveSpectrum){ //Background subtraction for IP analysis TH1D *incElecCent[kCentrality-1]; TH1D *charmCent[kCentrality-1]; TH1D *convCent[kCentrality-1]; TH1D *nonHFECent[kCentrality-1]; TH1D *subtractedCent[kCentrality-1]; TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr); CorrectFromTheWidth(measuredTH1Draw); if(fBeamType==1){ THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid(); for(Int_t icent = 1; icent < kCentrality-1; icent++){ sparseIncElec->GetAxis(0)->SetRange(icent,icent); incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr); CorrectFromTheWidth(incElecCent[icent-1]); } } TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400); rawspectra->cd(); rawspectra->SetLogy(); gStyle->SetOptStat(0); TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85); measuredTH1Draw->SetMarkerStyle(20); measuredTH1Draw->Draw(); measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9); lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum"); TH1D* htemp; Int_t* bins=new Int_t[2]; if(fIPanaHadronBgSubtract){ // Hadron background printf("Hadron background for IP analysis subtracted!\n"); if(fBeamType==0) { htemp = (TH1D *) fHadronEffbyIPcut->Projection(0); bins[0]=htemp->GetNbinsX(); } if(fBeamType==1) { htemp = (TH1D *) fHadronEffbyIPcut->Projection(0); bins[0]=htemp->GetNbinsX(); htemp = (TH1D *) fHadronEffbyIPcut->Projection(1); bins[1]=htemp->GetNbinsX(); } AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins); hbgContainer->SetGrid(fHadronEffbyIPcut); backgroundGrid->Multiply(hbgContainer,1); // draw raw hadron bg spectra TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr); CorrectFromTheWidth(hadronbg); hadronbg->SetMarkerColor(7); hadronbg->SetMarkerStyle(20); rawspectra->cd(); hadronbg->Draw("samep"); lRaw->AddEntry(hadronbg,"hadrons"); // subtract hadron contamination spectrumSubtracted->Add(backgroundGrid,-1.0); } if(fIPanaCharmBgSubtract){ // Charm background printf("Charm background for IP analysis subtracted!\n"); AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground(); // draw charm bg spectra TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr); CorrectFromTheWidth(charmbg); charmbg->SetMarkerColor(3); charmbg->SetMarkerStyle(20); rawspectra->cd(); charmbg->Draw("samep"); lRaw->AddEntry(charmbg,"charm elecs"); // subtract charm background spectrumSubtracted->Add(charmbgContainer,-1.0); if(fBeamType==1){ THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid(); for(Int_t icent = 1; icent < kCentrality-1; icent++){ sparseCharmElec->GetAxis(0)->SetRange(icent,icent); charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr); CorrectFromTheWidth(charmCent[icent-1]); } } } if(fIPanaConversionBgSubtract){ // Conversion background AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground(); // draw conversion bg spectra TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr); CorrectFromTheWidth(conversionbg); conversionbg->SetMarkerColor(4); conversionbg->SetMarkerStyle(20); rawspectra->cd(); conversionbg->Draw("samep"); lRaw->AddEntry(conversionbg,"conversion elecs"); // subtract conversion background spectrumSubtracted->Add(conversionbgContainer,-1.0); if(fBeamType==1){ THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid(); for(Int_t icent = 1; icent < kCentrality-1; icent++){ sparseconvElec->GetAxis(0)->SetRange(icent,icent); convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr); CorrectFromTheWidth(convCent[icent-1]); } } } if(fIPanaNonHFEBgSubtract){ // NonHFE background AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(); // draw Dalitz/dielectron bg spectra TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr); CorrectFromTheWidth(nonhfebg); nonhfebg->SetMarkerColor(6); nonhfebg->SetMarkerStyle(20); rawspectra->cd(); nonhfebg->Draw("samep"); lRaw->AddEntry(nonhfebg,"non-HF elecs"); // subtract Dalitz/dielectron background spectrumSubtracted->Add(nonHFEbgContainer,-1.0); if(fBeamType==1){ THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid(); for(Int_t icent = 1; icent < kCentrality-1; icent++){ sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent); nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr); CorrectFromTheWidth(nonHFECent[icent-1]); } } } TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr); CorrectFromTheWidth(rawbgsubtracted); rawbgsubtracted->SetMarkerStyle(24); rawspectra->cd(); lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum"); rawbgsubtracted->Draw("samep"); lRaw->Draw("SAME"); gPad->SetGrid(); //rawspectra->SaveAs("rawspectra.eps"); if(fBeamType==1){ THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid(); for(Int_t icent = 1; icent < kCentrality-1; icent++){ sparseSubtracted->GetAxis(0)->SetRange(icent,icent); subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr); CorrectFromTheWidth(subtractedCent[icent-1]); } TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85); TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800); centRaw->Divide(3,3); for(Int_t icent = 1; icent < kCentrality-1; icent++){ centRaw->cd(icent); gPad->SetLogx(); gPad->SetLogy(); incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.); incElecCent[icent-1]->Draw("p"); incElecCent[icent-1]->SetMarkerColor(1); incElecCent[icent-1]->SetMarkerStyle(20); charmCent[icent-1]->Draw("samep"); charmCent[icent-1]->SetMarkerColor(3); charmCent[icent-1]->SetMarkerStyle(20); convCent[icent-1]->Draw("samep"); convCent[icent-1]->SetMarkerColor(4); convCent[icent-1]->SetMarkerStyle(20); nonHFECent[icent-1]->Draw("samep"); nonHFECent[icent-1]->SetMarkerColor(6); nonHFECent[icent-1]->SetMarkerStyle(20); subtractedCent[icent-1]->Draw("samep"); subtractedCent[icent-1]->SetMarkerStyle(24); if(icent == 1){ lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum"); lCentRaw->AddEntry(charmCent[0],"charm elecs"); lCentRaw->AddEntry(convCent[0],"conversion elecs"); lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs"); lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum"); lCentRaw->Draw("SAME"); } } } delete[] bins; } else{ // Subtract spectrumSubtracted->Add(backgroundGrid,-1.0); } if(setBackground){ if(fBackground) delete fBackground; fBackground = backgroundGrid; } else delete backgroundGrid; if(fDebugLevel > 0) { Int_t ptprd; if(fBeamType==0) ptprd=0; if(fBeamType==1) ptprd=1; TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700); cbackgroundsubtraction->Divide(3,1); cbackgroundsubtraction->cd(1); gPad->SetLogy(); TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd); TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd); CorrectFromTheWidth(measuredTH1Daftersubstraction); CorrectFromTheWidth(measuredTH1Dbeforesubstraction); measuredTH1Daftersubstraction->SetStats(0); measuredTH1Daftersubstraction->SetTitle(""); measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]"); measuredTH1Daftersubstraction->SetMarkerStyle(25); measuredTH1Daftersubstraction->SetMarkerColor(kBlack); measuredTH1Daftersubstraction->SetLineColor(kBlack); measuredTH1Dbeforesubstraction->SetStats(0); measuredTH1Dbeforesubstraction->SetTitle(""); measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]"); measuredTH1Dbeforesubstraction->SetMarkerStyle(24); measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue); measuredTH1Dbeforesubstraction->SetLineColor(kBlue); measuredTH1Daftersubstraction->Draw(); measuredTH1Dbeforesubstraction->Draw("same"); TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89); legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p"); legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p"); legsubstraction->Draw("same"); cbackgroundsubtraction->cd(2); gPad->SetLogy(); TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone(); ratiomeasuredcontamination->SetName("ratiomeasuredcontamination"); ratiomeasuredcontamination->SetTitle(""); ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination"); ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratiomeasuredcontamination->Sumw2(); ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0); ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction); ratiomeasuredcontamination->SetStats(0); ratiomeasuredcontamination->SetMarkerStyle(26); ratiomeasuredcontamination->SetMarkerColor(kBlack); ratiomeasuredcontamination->SetLineColor(kBlack); for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){ ratiomeasuredcontamination->SetBinError(k+1,0.0); } ratiomeasuredcontamination->Draw("P"); cbackgroundsubtraction->cd(3); TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd); CorrectFromTheWidth(measuredTH1background); measuredTH1background->SetStats(0); measuredTH1background->SetTitle(""); measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]"); measuredTH1background->SetMarkerStyle(26); measuredTH1background->SetMarkerColor(kBlack); measuredTH1background->SetLineColor(kBlack); measuredTH1background->Draw(); if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps"); if(fBeamType==1) { TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700); cbackgrounde->Divide(2,1); TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89); TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89); THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid(); TAxis *cenaxisa = sparsesubtracted->GetAxis(0); THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid(); TAxis *cenaxisb = sparsebefore->GetAxis(0); Int_t nbbin = cenaxisb->GetNbins(); Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29}; Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20}; for(Int_t binc = 0; binc < nbbin; binc++){ TString titlee("BackgroundSubtraction_centrality_bin_"); titlee += binc; TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700); cbackground->Divide(2,1); cbackground->cd(1); gPad->SetLogy(); cenaxisa->SetRange(binc+1,binc+1); cenaxisb->SetRange(binc+1,binc+1); TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1); TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1); CorrectFromTheWidth(aftersubstraction); CorrectFromTheWidth(beforesubstraction); aftersubstraction->SetStats(0); aftersubstraction->SetTitle((const char*)titlee); aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]"); aftersubstraction->SetMarkerStyle(25); aftersubstraction->SetMarkerColor(kBlack); aftersubstraction->SetLineColor(kBlack); beforesubstraction->SetStats(0); beforesubstraction->SetTitle((const char*)titlee); beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]"); beforesubstraction->SetMarkerStyle(24); beforesubstraction->SetMarkerColor(kBlue); beforesubstraction->SetLineColor(kBlue); aftersubstraction->Draw(); beforesubstraction->Draw("same"); TLegend *lega = new TLegend(0.4,0.6,0.89,0.89); lega->AddEntry(beforesubstraction,"With hadron contamination","p"); lega->AddEntry(aftersubstraction,"Without hadron contamination ","p"); lega->Draw("same"); cbackgrounde->cd(1); gPad->SetLogy(); TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone(); aftersubtractionn->SetMarkerStyle(stylee[binc]); aftersubtractionn->SetMarkerColor(colorr[binc]); if(binc==0) aftersubtractionn->Draw(); else aftersubtractionn->Draw("same"); legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p"); cbackgrounde->cd(2); gPad->SetLogy(); TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone(); aftersubtractionng->SetMarkerStyle(stylee[binc]); aftersubtractionng->SetMarkerColor(colorr[binc]); if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]); if(binc==0) aftersubtractionng->Draw(); else aftersubtractionng->Draw("same"); legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p"); cbackground->cd(2); TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone(); ratiocontamination->SetName("ratiocontamination"); ratiocontamination->SetTitle((const char*)titlee); ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination"); ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratiocontamination->Add(aftersubstraction,-1.0); ratiocontamination->Divide(beforesubstraction); Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins(); for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) { ratiocontamination->SetBinError(nbinpt+1,0.0); } ratiocontamination->SetStats(0); ratiocontamination->SetMarkerStyle(26); ratiocontamination->SetMarkerColor(kBlack); ratiocontamination->SetLineColor(kBlack); ratiocontamination->Draw("P"); } cbackgrounde->cd(1); legtotal->Draw("same"); cbackgrounde->cd(2); legtotalg->Draw("same"); cenaxisa->SetRange(0,nbbin); cenaxisb->SetRange(0,nbbin); if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps"); } } return spectrumSubtracted; } //____________________________________________________________ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){ // // calculate charm background // Int_t ptpr = 0; Int_t nDim = 1; if(fBeamType==0) { ptpr=0; } if(fBeamType==1) { ptpr=1; nDim=2; } Double_t evtnorm=0; if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]); AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC); if(!mcContainer){ AliError("MC Container not available"); return NULL; } if(!fCorrelation){ AliError("No Correlation map available"); return NULL; } AliCFDataGrid *charmBackgroundGrid= 0x0; charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0); Int_t* bins=new Int_t[2]; bins[0]=charmbgaftertofpid->GetNbinsX(); if(fBeamType==1) { bins[0]=12; charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1); bins[1]=charmbgaftertofpid->GetNbinsX(); } AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins); ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr); charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.); Int_t *nBinpp=new Int_t[1]; Int_t* binspp=new Int_t[1]; binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins Int_t *nBinPbPb=new Int_t[2]; Int_t* binsPbPb=new Int_t[2]; binsPbPb[1]=charmbgaftertofpid->GetNbinsX(); // number of pt bins binsPbPb[0]=12; Int_t looppt=binspp[0]; if(fBeamType==1) looppt=binsPbPb[1]; for(Long_t iBin=1; iBin<= looppt;iBin++){ if(fBeamType==0) { nBinpp[0]=iBin; charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm); charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm); } if(fBeamType==1) { // loop over centrality for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){ nBinPbPb[0]=iiBin; nBinPbPb[1]=iBin; Double_t evtnormPbPb=0; if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]); charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb); charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb); } } } TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr); AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins); weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr); charmBackgroundGrid->Multiply(weightedCharmContainer,1.); TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr); // Efficiency (set efficiency to 1 for only folding) AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(0,0); // Folding if(fBeamType==0)nDim = 1; if(fBeamType==1)nDim = 2; AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid()); folding.SetMaxNumberOfIterations(1); folding.Unfold(); // Results THnSparse* result1= folding.GetEstMeasured(); // folded spectra THnSparse* result=(THnSparse*)result1->Clone(); charmBackgroundGrid->SetGrid(result); TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr); //Charm background evaluation plots TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600); cCharmBgEval->Divide(3,1); cCharmBgEval->cd(1); if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm); if(fBeamType==1) { Double_t evtnormPbPb=0; for(Int_t kCentr=0;kCentrScale(evtnormPbPb); } CorrectFromTheWidth(charmbgaftertofpid); charmbgaftertofpid->SetMarkerStyle(25); charmbgaftertofpid->Draw("p"); charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events"); charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)"); gPad->SetLogy(); CorrectFromTheWidth(charmbgafteripcut); charmbgafteripcut->SetMarkerStyle(24); charmbgafteripcut->Draw("samep"); CorrectFromTheWidth(charmbgafterweight); charmbgafterweight->SetMarkerStyle(24); charmbgafterweight->SetMarkerColor(4); charmbgafterweight->Draw("samep"); CorrectFromTheWidth(charmbgafterfolding); charmbgafterfolding->SetMarkerStyle(24); charmbgafterfolding->SetMarkerColor(2); charmbgafterfolding->Draw("samep"); cCharmBgEval->cd(2); parametrizedcharmpidipeff->SetMarkerStyle(24); parametrizedcharmpidipeff->Draw("p"); parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)"); cCharmBgEval->cd(3); charmweightingfc->SetMarkerStyle(24); charmweightingfc->Draw("p"); charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron"); charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)"); cCharmBgEval->cd(1); TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89); legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p"); legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p"); legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p"); legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p"); legcharmbg->Draw("same"); cCharmBgEval->cd(2); TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89); legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p"); legcharmbg2->Draw("same"); CorrectStatErr(charmBackgroundGrid); if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps"); delete[] bins; delete[] nBinpp; delete[] binspp; delete[] nBinPbPb; delete[] binsPbPb; if(fUnfoldBG) UnfoldBG(charmBackgroundGrid); return charmBackgroundGrid; } //____________________________________________________________ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){ // // calculate conversion background // Double_t evtnorm[1] = {0.0}; if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); printf("check event!!! %lf \n",evtnorm[0]); AliCFContainer *backgroundContainer = 0x0; if(fNonHFEsyst){ backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone(); for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent } } else{ // Background Estimate backgroundContainer = GetContainer(kMCWeightedContainerConversionESD); } if(!backgroundContainer){ AliError("MC background container not found"); return NULL; } Int_t stepbackground = 3; AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground); Int_t *nBinpp=new Int_t[1]; Int_t* binspp=new Int_t[1]; binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins Int_t *nBinPbPb=new Int_t[2]; Int_t* binsPbPb=new Int_t[2]; binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins binsPbPb[0]=12; Int_t looppt=binspp[0]; if(fBeamType==1) looppt=binsPbPb[1]; for(Long_t iBin=1; iBin<= looppt;iBin++){ if(fBeamType==0) { nBinpp[0]=iBin; backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]); backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]); } if(fBeamType==1) { // loop over centrality for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){ nBinPbPb[0]=iiBin; nBinPbPb[1]=iBin; Double_t evtnormPbPb=0; if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]); backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb); backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb); } } } //end of workaround for statistical errors AliCFDataGrid *weightedConversionContainer; if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp); else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb); weightedConversionContainer->SetGrid(GetPIDxIPEff(2)); backgroundGrid->Multiply(weightedConversionContainer,1.0); delete[] nBinpp; delete[] binspp; delete[] nBinPbPb; delete[] binsPbPb; return backgroundGrid; } //____________________________________________________________ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){ // // calculate non-HFE background // Double_t evtnorm[1] = {0.0}; if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); printf("check event!!! %lf \n",evtnorm[0]); AliCFContainer *backgroundContainer = 0x0; if(fNonHFEsyst){ backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone(); for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ if(iSource == 1) backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA else backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]); } } else{ // Background Estimate backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD); } if(!backgroundContainer){ AliError("MC background container not found"); return NULL; } Int_t stepbackground = 3; AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground); Int_t *nBinpp=new Int_t[1]; Int_t* binspp=new Int_t[1]; binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins Int_t *nBinPbPb=new Int_t[2]; Int_t* binsPbPb=new Int_t[2]; binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins binsPbPb[0]=12; Int_t looppt=binspp[0]; if(fBeamType==1) looppt=binsPbPb[1]; for(Long_t iBin=1; iBin<= looppt;iBin++){ if(fBeamType==0) { nBinpp[0]=iBin; backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]); backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]); } if(fBeamType==1) { for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){ nBinPbPb[0]=iiBin; nBinPbPb[1]=iBin; Double_t evtnormPbPb=0; if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]); backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb); backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb); } } } //end of workaround for statistical errors AliCFDataGrid *weightedNonHFEContainer; if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp); else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb); weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3)); backgroundGrid->Multiply(weightedNonHFEContainer,1.0); delete[] nBinpp; delete[] binspp; delete[] nBinPbPb; delete[] binsPbPb; return backgroundGrid; } //____________________________________________________________ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){ // // Apply TPC pid efficiency correction from parametrisation // // Data in the right format AliCFDataGrid *dataGrid = 0x0; if(bgsubpectrum) { dataGrid = bgsubpectrum; } else { AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData); } AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone(); result->SetName("ParametrizedEfficiencyBefore"); THnSparse *h = result->GetGrid(); Int_t nbdimensions = h->GetNdimensions(); //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions); AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone(); dataContainerbis->Add(dataContainerbis,-1.0); Int_t* coord = new Int_t[nbdimensions]; memset(coord, 0, sizeof(Int_t) * nbdimensions); Double_t* points = new Double_t[nbdimensions]; ULong64_t nEntries = h->GetNbins(); for (ULong64_t i = 0; i < nEntries; ++i) { Double_t value = h->GetBinContent(i, coord); //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData); //printf("Value %f, and valuecontainer %f\n",value,valuecontainer); // Get the bin co-ordinates given an coord for (Int_t j = 0; j < nbdimensions; ++j) points[j] = h->GetAxis(j)->GetBinCenter(coord[j]); if (!fEfficiencyFunction->IsInside(points)) continue; TF1::RejectPoint(kFALSE); // Evaulate function at points Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL); //printf("Value efficiency is %f\n",valueEfficiency); if(valueEfficiency > 0.0) { h->SetBinContent(coord,value/valueEfficiency); dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency); } Double_t error = h->GetBinError(i); h->SetBinError(coord,error/valueEfficiency); dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency); } delete[] coord; delete[] points; AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData); if(fDebugLevel > 0) { TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700); cEfficiencyParametrized->Divide(2,1); cEfficiencyParametrized->cd(1); TH1D *afterE = (TH1D *) resultt->Project(0); TH1D *beforeE = (TH1D *) dataGrid->Project(0); CorrectFromTheWidth(afterE); CorrectFromTheWidth(beforeE); afterE->SetStats(0); afterE->SetTitle(""); afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]"); afterE->SetMarkerStyle(25); afterE->SetMarkerColor(kBlack); afterE->SetLineColor(kBlack); beforeE->SetStats(0); beforeE->SetTitle(""); beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]"); beforeE->SetMarkerStyle(24); beforeE->SetMarkerColor(kBlue); beforeE->SetLineColor(kBlue); gPad->SetLogy(); afterE->Draw(); beforeE->Draw("same"); TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89); legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p"); legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p"); legefficiencyparametrized->Draw("same"); cEfficiencyParametrized->cd(2); fEfficiencyFunction->Draw(); //cEfficiencyParametrized->cd(3); //TH1D *ratioefficiency = (TH1D *) beforeE->Clone(); //ratioefficiency->Divide(afterE); //ratioefficiency->Draw(); if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps"); } return resultt; } //____________________________________________________________ AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){ // // Apply TPC pid efficiency correction from V0 // AliCFContainer *v0Container = GetContainer(kDataContainerV0); if(!v0Container){ AliError("V0 Container not available"); return NULL; } // Efficiency AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container); efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0); // Data in the right format AliCFDataGrid *dataGrid = 0x0; if(bgsubpectrum) { dataGrid = bgsubpectrum; } else { AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData); } // Correct AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone(); result->ApplyEffCorrection(*efficiencyD); if(fDebugLevel > 0) { Int_t ptpr; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700); cV0Efficiency->Divide(2,1); cV0Efficiency->cd(1); TH1D *afterE = (TH1D *) result->Project(ptpr); TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr); afterE->SetStats(0); afterE->SetTitle(""); afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]"); afterE->SetMarkerStyle(25); afterE->SetMarkerColor(kBlack); afterE->SetLineColor(kBlack); beforeE->SetStats(0); beforeE->SetTitle(""); beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]"); beforeE->SetMarkerStyle(24); beforeE->SetMarkerColor(kBlue); beforeE->SetLineColor(kBlue); afterE->Draw(); beforeE->Draw("same"); TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89); legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p"); legV0efficiency->AddEntry(afterE,"After Efficiency correction","p"); legV0efficiency->Draw("same"); cV0Efficiency->cd(2); TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr); efficiencyDproj->SetTitle(""); efficiencyDproj->SetStats(0); efficiencyDproj->SetMarkerStyle(25); efficiencyDproj->Draw(); if(fBeamType==1) { TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700); cV0Efficiencye->Divide(2,1); TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89); TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89); THnSparseF* sparseafter = (THnSparseF *) result->GetGrid(); TAxis *cenaxisa = sparseafter->GetAxis(0); THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid(); TAxis *cenaxisb = sparsebefore->GetAxis(0); THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid(); TAxis *cenaxisc = efficiencya->GetAxis(0); Int_t nbbin = cenaxisb->GetNbins(); Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29}; Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20}; for(Int_t binc = 0; binc < nbbin; binc++){ TString titlee("V0Efficiency_centrality_bin_"); titlee += binc; TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700); ccV0Efficiency->Divide(2,1); ccV0Efficiency->cd(1); gPad->SetLogy(); cenaxisa->SetRange(binc+1,binc+1); cenaxisb->SetRange(binc+1,binc+1); cenaxisc->SetRange(binc+1,binc+1); TH1D *aftere = (TH1D *) sparseafter->Projection(1); TH1D *beforee = (TH1D *) sparsebefore->Projection(1); CorrectFromTheWidth(aftere); CorrectFromTheWidth(beforee); aftere->SetStats(0); aftere->SetTitle((const char*)titlee); aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]"); aftere->SetMarkerStyle(25); aftere->SetMarkerColor(kBlack); aftere->SetLineColor(kBlack); beforee->SetStats(0); beforee->SetTitle((const char*)titlee); beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]"); beforee->SetMarkerStyle(24); beforee->SetMarkerColor(kBlue); beforee->SetLineColor(kBlue); aftere->Draw(); beforee->Draw("same"); TLegend *lega = new TLegend(0.4,0.6,0.89,0.89); lega->AddEntry(beforee,"Before correction","p"); lega->AddEntry(aftere,"After correction","p"); lega->Draw("same"); cV0Efficiencye->cd(1); gPad->SetLogy(); TH1D *afteree = (TH1D *) aftere->Clone(); afteree->SetMarkerStyle(stylee[binc]); afteree->SetMarkerColor(colorr[binc]); if(binc==0) afteree->Draw(); else afteree->Draw("same"); legtotal->AddEntry(afteree,(const char*) titlee,"p"); cV0Efficiencye->cd(2); gPad->SetLogy(); TH1D *aftereeu = (TH1D *) aftere->Clone(); aftereeu->SetMarkerStyle(stylee[binc]); aftereeu->SetMarkerColor(colorr[binc]); if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]); if(binc==0) aftereeu->Draw(); else aftereeu->Draw("same"); legtotalg->AddEntry(aftereeu,(const char*) titlee,"p"); ccV0Efficiency->cd(2); TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1); efficiencyDDproj->SetTitle(""); efficiencyDDproj->SetStats(0); efficiencyDDproj->SetMarkerStyle(25); efficiencyDDproj->Draw(); } cV0Efficiencye->cd(1); legtotal->Draw("same"); cV0Efficiencye->cd(2); legtotalg->Draw("same"); cenaxisa->SetRange(0,nbbin); cenaxisb->SetRange(0,nbbin); cenaxisc->SetRange(0,nbbin); if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps"); } } return result; } //____________________________________________________________ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){ // // Unfold and eventually correct for efficiency the bgsubspectrum // AliCFContainer *mcContainer = GetContainer(kMCContainerMC); if(!mcContainer){ AliError("MC Container not available"); return NULL; } if(!fCorrelation){ AliError("No Correlation map available"); return NULL; } // Data AliCFDataGrid *dataGrid = 0x0; if(bgsubpectrum) { dataGrid = bgsubpectrum; } else { AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData); } // Guessed AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding); THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid(); // Efficiency AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(fStepMC,fStepTrue); if(!fBeauty2ndMethod) { // Consider parameterized IP cut efficiency if(!fInclusiveSpectrum){ Int_t* bins=new Int_t[1]; bins[0]=35; AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); beffContainer->SetGrid(GetBeautyIPEff(kTRUE)); efficiencyD->Multiply(beffContainer,1); } } // Unfold AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse); if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors(); unfolding.SetMaxNumberOfIterations(fNumberOfIterations); if(fNRandomIter > 0) unfolding.SetNRandomIterations(fNRandomIter); if(fSetSmoothing) unfolding.UseSmoothing(); unfolding.Unfold(); // Results THnSparse* result = unfolding.GetUnfolded(); THnSparse* residual = unfolding.GetEstMeasured(); TList *listofresults = new TList; listofresults->SetOwner(); listofresults->AddAt((THnSparse*)result->Clone(),0); listofresults->AddAt((THnSparse*)residual->Clone(),1); if(fDebugLevel > 0) { Int_t ptpr; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TCanvas * cresidual = new TCanvas("residual","residual",1000,700); cresidual->Divide(2,1); cresidual->cd(1); gPad->SetLogy(); TGraphErrors* residualspectrumD = Normalize(residual); if(!residualspectrumD) { AliError("Number of Events not set for the normalization"); return NULL; } residualspectrumD->SetTitle(""); residualspectrumD->GetYaxis()->SetTitleOffset(1.5); residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); residualspectrumD->SetMarkerStyle(26); residualspectrumD->SetMarkerColor(kBlue); residualspectrumD->SetLineColor(kBlue); residualspectrumD->Draw("AP"); AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone()); dataGridBis->SetName("dataGridBis"); TGraphErrors* measuredspectrumD = Normalize(dataGridBis); measuredspectrumD->SetTitle(""); measuredspectrumD->GetYaxis()->SetTitleOffset(1.5); measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0); measuredspectrumD->SetMarkerStyle(25); measuredspectrumD->SetMarkerColor(kBlack); measuredspectrumD->SetLineColor(kBlack); measuredspectrumD->Draw("P"); TLegend *legres = new TLegend(0.4,0.6,0.89,0.89); legres->AddEntry(residualspectrumD,"Residual","p"); legres->AddEntry(measuredspectrumD,"Measured","p"); legres->Draw("same"); cresidual->cd(2); TH1D *residualTH1D = residual->Projection(ptpr); TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr); TH1D* ratioresidual = (TH1D*)residualTH1D->Clone(); ratioresidual->SetName("ratioresidual"); ratioresidual->SetTitle(""); ratioresidual->GetYaxis()->SetTitle("Folded/Measured"); ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]"); ratioresidual->Divide(residualTH1D,measuredTH1D,1,1); ratioresidual->SetStats(0); ratioresidual->Draw(); if(fWriteToFile) cresidual->SaveAs("Unfolding.eps"); } return listofresults; } //____________________________________________________________ AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){ // // Apply unfolding and efficiency correction together to bgsubspectrum // AliCFContainer *mcContainer = GetContainer(kMCContainerESD); if(!mcContainer){ AliError("MC Container not available"); return NULL; } // Efficiency AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(fStepMC,fStepTrue); if(!fBeauty2ndMethod) { // Consider parameterized IP cut efficiency if(!fInclusiveSpectrum){ Int_t* bins=new Int_t[1]; bins[0]=35; AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); beffContainer->SetGrid(GetBeautyIPEff(kFALSE)); efficiencyD->Multiply(beffContainer,1); } } // Data in the right format AliCFDataGrid *dataGrid = 0x0; if(bgsubpectrum) { dataGrid = bgsubpectrum; } else { AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); return NULL; } dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData); } // Correct AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone(); result->ApplyEffCorrection(*efficiencyD); if(fDebugLevel > 0) { Int_t ptpr; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; printf("Step MC: %d\n",fStepMC); printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack); printf("Step MC true: %d\n",fStepTrue); AliCFEffGrid *efficiencymcPID = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack); AliCFEffGrid *efficiencymctrackinggeo = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue); AliCFEffGrid *efficiencymcall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue); AliCFEffGrid *efficiencyesdall = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue); TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700); cefficiency->cd(1); TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr); efficiencymcPIDD->SetTitle(""); efficiencymcPIDD->SetStats(0); efficiencymcPIDD->SetMarkerStyle(25); efficiencymcPIDD->Draw(); TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr); efficiencymctrackinggeoD->SetTitle(""); efficiencymctrackinggeoD->SetStats(0); efficiencymctrackinggeoD->SetMarkerStyle(26); efficiencymctrackinggeoD->Draw("same"); TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr); efficiencymcallD->SetTitle(""); efficiencymcallD->SetStats(0); efficiencymcallD->SetMarkerStyle(27); efficiencymcallD->Draw("same"); TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr); efficiencyesdallD->SetTitle(""); efficiencyesdallD->SetStats(0); efficiencyesdallD->SetMarkerStyle(24); efficiencyesdallD->Draw("same"); TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89); legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p"); legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p"); legeff->AddEntry(efficiencymcallD,"Overall efficiency","p"); legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p"); legeff->Draw("same"); if(fBeamType==1) { THnSparseF* sparseafter = (THnSparseF *) result->GetGrid(); TAxis *cenaxisa = sparseafter->GetAxis(0); THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid(); TAxis *cenaxisb = sparsebefore->GetAxis(0); THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid(); TAxis *cenaxisc = efficiencya->GetAxis(0); //Int_t nbbin = cenaxisb->GetNbins(); //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29}; //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20}; for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){ TString titlee("Efficiency_centrality_bin_"); titlee += fLowBoundaryCentralityBinAtTheEnd[binc]; titlee += "_"; titlee += fHighBoundaryCentralityBinAtTheEnd[binc]; TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700); cefficiencye->Divide(2,1); cefficiencye->cd(1); gPad->SetLogy(); cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]); cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]); TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr); TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr); CorrectFromTheWidth(afterefficiency); CorrectFromTheWidth(beforeefficiency); afterefficiency->SetStats(0); afterefficiency->SetTitle((const char*)titlee); afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]"); afterefficiency->SetMarkerStyle(25); afterefficiency->SetMarkerColor(kBlack); afterefficiency->SetLineColor(kBlack); beforeefficiency->SetStats(0); beforeefficiency->SetTitle((const char*)titlee); beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]"); beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]"); beforeefficiency->SetMarkerStyle(24); beforeefficiency->SetMarkerColor(kBlue); beforeefficiency->SetLineColor(kBlue); afterefficiency->Draw(); beforeefficiency->Draw("same"); TLegend *lega = new TLegend(0.4,0.6,0.89,0.89); lega->AddEntry(beforeefficiency,"Before efficiency correction","p"); lega->AddEntry(afterefficiency,"After efficiency correction","p"); lega->Draw("same"); cefficiencye->cd(2); cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1); TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr); efficiencyDDproj->SetTitle(""); efficiencyDDproj->SetStats(0); efficiencyDDproj->SetMarkerStyle(25); efficiencyDDproj->SetMarkerColor(2); efficiencyDDproj->Draw(); cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]); TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr); efficiencyDDproja->SetTitle(""); efficiencyDDproja->SetStats(0); efficiencyDDproja->SetMarkerStyle(26); efficiencyDDproja->SetMarkerColor(4); efficiencyDDproja->Draw("same"); } } if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps"); } return result; } //__________________________________________________________________________________ TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const { // // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2} // Give the final pt spectrum to be compared // if(fNEvents[i] > 0) { Int_t ptpr = 0; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TH1D* projection = spectrum->Projection(ptpr); CorrectFromTheWidth(projection); TGraphErrors *graphError = NormalizeTH1(projection,i); return graphError; } return 0x0; } //__________________________________________________________________________________ TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const { // // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2} // Give the final pt spectrum to be compared // if(fNEvents[i] > 0) { Int_t ptpr=0; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; TH1D* projection = (TH1D *) spectrum->Project(ptpr); CorrectFromTheWidth(projection); TGraphErrors *graphError = NormalizeTH1(projection,i); return graphError; } return 0x0; } //__________________________________________________________________________________ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const { // // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2} // Give the final pt spectrum to be compared // Double_t chargecoefficient = 0.5; if(fChargeChoosen != 0) chargecoefficient = 1.0; Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6; printf("Normalizing Eta Range %f\n", etarange); if(fNEvents[i] > 0) { TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX()); Double_t p = 0, dp = 0; Int_t point = 1; Double_t n = 0, dN = 0; Double_t nCorr = 0, dNcorr = 0; //Double_t errdN = 0, errdp = 0; Double_t errdN = 0; for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){ point = ibin - input->GetXaxis()->GetFirst(); p = input->GetXaxis()->GetBinCenter(ibin); dp = input->GetXaxis()->GetBinWidth(ibin)/2.; n = input->GetBinContent(ibin); dN = input->GetBinError(ibin); // New point nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n; errdN = 1./(2. * TMath::Pi() * p); //errdp = 1./(2. * TMath::Pi() * p*p) * n; //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp); dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN); spectrumNormalized->SetPoint(point, p, nCorr); spectrumNormalized->SetPointError(point, dp, dNcorr); } spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]"); spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}"); spectrumNormalized->SetMarkerStyle(22); spectrumNormalized->SetMarkerColor(kBlue); spectrumNormalized->SetLineColor(kBlue); return spectrumNormalized; } return 0x0; } //__________________________________________________________________________________ TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const { // // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2} // Give the final pt spectrum to be compared // Double_t chargecoefficient = 0.5; if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0; Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6; printf("Normalizing Eta Range %f\n", etarange); if(normalization > 0) { TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX()); Double_t p = 0, dp = 0; Int_t point = 1; Double_t n = 0, dN = 0; Double_t nCorr = 0, dNcorr = 0; //Double_t errdN = 0, errdp = 0; Double_t errdN = 0; for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){ point = ibin - input->GetXaxis()->GetFirst(); p = input->GetXaxis()->GetBinCenter(ibin); //dp = input->GetXaxis()->GetBinWidth(ibin)/2.; n = input->GetBinContent(ibin); dN = input->GetBinError(ibin); // New point nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n; errdN = 1./(2. * TMath::Pi() * p); //errdp = 1./(2. * TMath::Pi() * p*p) * n; //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp); dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN); spectrumNormalized->SetPoint(point, p, nCorr); spectrumNormalized->SetPointError(point, dp, dNcorr); } spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]"); spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}"); spectrumNormalized->SetMarkerStyle(22); spectrumNormalized->SetMarkerColor(kBlue); spectrumNormalized->SetLineColor(kBlue); return spectrumNormalized; } return 0x0; } //____________________________________________________________ void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){ // // Set the container for a given step to the // if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1); fCFContainers->AddAt(cont, type); } //____________________________________________________________ AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){ // // Get Correction Framework Container for given type // if(!fCFContainers) return NULL; return dynamic_cast(fCFContainers->At(type)); } //____________________________________________________________ AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) { // // Slice bin for a given source of electron // nDim is the number of dimension the corrections are done // dimensions are the definition of the dimensions // source is if we want to keep only one MC source (-1 means we don't cut on the MC source) // positivenegative if we want to keep positive (1) or negative (0) or both (-1) // centrality (-1 means we do not cut on centrality) // Double_t *varMin = new Double_t[container->GetNVar()], *varMax = new Double_t[container->GetNVar()]; Double_t *binLimits; for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){ binLimits = new Double_t[container->GetNBins(ivar)+1]; container->GetBinLimits(ivar,binLimits); varMin[ivar] = binLimits[0]; varMax[ivar] = binLimits[container->GetNBins(ivar)]; // source if(ivar == 4){ if((source>= 0) && (sourceGetNBins(ivar))) { varMin[ivar] = container->GetAxis(4,0)->GetBinLowEdge(container->GetAxis(4,0)->FindBin(binLimits[source])); varMax[ivar] = container->GetAxis(4,0)->GetBinUpEdge(container->GetAxis(4,0)->FindBin(binLimits[source])); } } // charge if(ivar == 3) { if(charge != kAllCharge){ varMin[ivar] = container->GetAxis(3,0)->GetBinLowEdge(container->GetAxis(3,0)->FindBin(charge)); varMax[ivar] = container->GetAxis(3,0)->GetBinUpEdge(container->GetAxis(3,0)->FindBin(charge)); } } // eta if(ivar == 1){ for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic))); if(fEtaSelected){ varMin[ivar] = fEtaRange[0]; varMax[ivar] = fEtaRange[1]; } } if(fEtaSelected){ fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0])); fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1])); AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0])); } if(ivar == 5){ if((centralitylow>= 0) && (centralitylowGetNBins(ivar)) && (centralityhigh>= 0) && (centralityhighGetNBins(ivar))) { varMin[ivar] = binLimits[centralitylow]; varMax[ivar] = binLimits[centralityhigh]; TAxis *axistest = container->GetAxis(5,0); printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins()); printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]); Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow])); Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh])); printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality); } } delete[] binLimits; } AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax); AddTemporaryObject(k); delete[] varMin; delete[] varMax; return k; } //_________________________________________________________________________ THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const { // // Slice correlation // Int_t ndimensions = correlationmatrix->GetNdimensions(); //printf("Number of dimension %d correlation map\n",ndimensions); if(ndimensions < (2*nDim)) { AliError("Problem in the dimensions"); return NULL; } // Cut in centrality is centrality > -1 if((centralitylow >=0) && (centralityhigh >=0)) { TAxis *axiscentrality0 = correlationmatrix->GetAxis(5); TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.))); Int_t bins0 = axiscentrality0->GetNbins(); Int_t bins1 = axiscentrality1->GetNbins(); printf("Number of centrality bins: %d and %d\n",bins0,bins1); if(bins0 != bins1) { AliError("Problem in the dimensions"); return NULL; } if((centralitylow>= 0) && (centralitylow= 0) && (centralityhighSetRangeUser(centralitylow,centralityhigh); axiscentrality1->SetRangeUser(centralitylow,centralityhigh); Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow)); Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh)); Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow)); Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh)); printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0); printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1); TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700); ctestcorrelation->cd(1); TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.))); projection->Draw("colz"); } } Int_t ndimensionsContainer = (Int_t) ndimensions/2; //printf("Number of dimension %d container\n",ndimensionsContainer); Int_t *dim = new Int_t[nDim*2]; for(Int_t iter=0; iter < nDim; iter++){ dim[iter] = dimensions[iter]; dim[iter+nDim] = ndimensionsContainer + dimensions[iter]; //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]); } THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim); delete[] dim; return k; } //___________________________________________________________________________ void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const { // // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1} // TAxis *axis = h1->GetXaxis(); Int_t nbinX = h1->GetNbinsX(); for(Int_t i = 1; i <= nbinX; i++) { Double_t width = axis->GetBinWidth(i); Double_t content = h1->GetBinContent(i); Double_t error = h1->GetBinError(i); h1->SetBinContent(i,content/width); h1->SetBinError(i,error/width); } } //___________________________________________________________________________ void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const { // // Correct statistical error // TH1D *h1 = (TH1D*)backgroundGrid->Project(0); Int_t nbinX = h1->GetNbinsX(); Int_t bins[1]; for(Long_t i = 1; i <= nbinX; i++) { bins[0] = i; Float_t content = h1->GetBinContent(i); if(content>0){ Float_t error = TMath::Sqrt(content); backgroundGrid->SetElementError(bins, error); } } } //____________________________________________________________ void AliHFEspectrum::AddTemporaryObject(TObject *o){ // // Emulate garbage collection on class level // if(!fTemporaryObjects) { fTemporaryObjects= new TList; fTemporaryObjects->SetOwner(); } fTemporaryObjects->Add(o); } //____________________________________________________________ void AliHFEspectrum::ClearObject(TObject *o){ // // Do a safe deletion // if(fTemporaryObjects){ if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o); delete o; } else{ // just delete delete o; } } //_________________________________________________________________________ TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) { AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step); return data; } //_________________________________________________________________________ TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){ // // Create efficiency grid and calculate efficiency // of step to step0 // TString name("eff"); name += step; name+= step0; AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c); eff->CalculateEfficiency(step,step0); return eff; } //____________________________________________________________________________ THnSparse* AliHFEspectrum::GetCharmWeights(){ // // Measured D->e based weighting factors // const Int_t nDimpp=1; Int_t nBinpp[nDimpp] = {35}; 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.}; const Int_t nDimPbPb=2; Int_t nBinPbPb[nDimPbPb] = {11,35}; Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; Int_t loopcentr=1; Int_t looppt=nBinpp[0]; if(fBeamType==0) { fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp); fWeightCharm->SetBinEdges(0, ptbinning1); } if(fBeamType==1) { fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb); fWeightCharm->SetBinEdges(1, ptbinning1); fWeightCharm->SetBinEdges(0, kCentralityRange); loopcentr=nBinPbPb[0]; } // Weighting factor for pp 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}; // Weighting factor for PbPb (0-20%) //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}; // Weighting factor for PbPb (40-80%) //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}; //points Double_t pt[1]; Double_t contents[2]; for(int icentr=0; icentrFill(contents,weight[i]); } } Int_t nDimSparse = fWeightCharm->GetNdimensions(); Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); nBins *= binsvar[iVar]; } Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) // loop that sets 0 error in each bin for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere } delete[] binsvar; delete[] binfill; return fWeightCharm; } //____________________________________________________________________________ void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){ // TOF PID efficiencies Int_t ptpr=0; if(fBeamType==0) ptpr=0; if(fBeamType==1) ptpr=1; Int_t loopcentr=1; const Int_t nCentralitybinning=11; //number of centrality bins if(fBeamType==1) { loopcentr=nCentralitybinning; } TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.); TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.); TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0); TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600); cefficiencyParamtof->cd(); AliCFContainer *mccontainermcD = 0x0; AliCFContainer *mccontaineresdD = 0x0; TH1D* efficiencysigTOFPIDD[nCentralitybinning]; TH1D* efficiencyTOFPIDD[nCentralitybinning]; TH1D* efficiencysigesdTOFPIDD[nCentralitybinning]; TH1D* efficiencyesdTOFPIDD[nCentralitybinning]; Int_t source = -1; //get parameterized TOF PID efficiencies for(int icentr=0; icentrCalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies // mb sample for double check if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1); AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD); efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies // mb sample with reconstructed variables if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1); AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD); efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies // mb sample with reconstructed variables if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1); AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD); efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies //fill histo efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr); efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr); efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr); efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr); efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr)); efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr)); efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr)); efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr)); //fit (mc enhenced sample) fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R"); efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency"); fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid"); //fit (esd enhenced sample) efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R"); efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency"); fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid"); } // draw (for PbPb, only 1st bin) //sig mc efficiencysigTOFPIDD[0]->SetTitle(""); efficiencysigTOFPIDD[0]->SetStats(0); efficiencysigTOFPIDD[0]->SetMarkerStyle(25); efficiencysigTOFPIDD[0]->SetMarkerColor(2); efficiencysigTOFPIDD[0]->SetLineColor(2); efficiencysigTOFPIDD[0]->Draw(); //mb mc efficiencyTOFPIDD[0]->SetTitle(""); efficiencyTOFPIDD[0]->SetStats(0); efficiencyTOFPIDD[0]->SetMarkerStyle(24); efficiencyTOFPIDD[0]->SetMarkerColor(4); efficiencyTOFPIDD[0]->SetLineColor(4); efficiencyTOFPIDD[0]->Draw("same"); //sig esd efficiencysigesdTOFPIDD[0]->SetTitle(""); efficiencysigesdTOFPIDD[0]->SetStats(0); efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25); efficiencysigesdTOFPIDD[0]->SetMarkerColor(3); efficiencysigesdTOFPIDD[0]->SetLineColor(3); efficiencysigesdTOFPIDD[0]->Draw("same"); //mb esd efficiencyesdTOFPIDD[0]->SetTitle(""); efficiencyesdTOFPIDD[0]->SetStats(0); efficiencyesdTOFPIDD[0]->SetMarkerStyle(25); efficiencyesdTOFPIDD[0]->SetMarkerColor(1); efficiencyesdTOFPIDD[0]->SetLineColor(1); efficiencyesdTOFPIDD[0]->Draw("same"); //signal mc fit if(fEfficiencyTOFPIDD[0]){ fEfficiencyTOFPIDD[0]->SetLineColor(2); fEfficiencyTOFPIDD[0]->Draw("same"); } //mb esd fit if(fEfficiencyesdTOFPIDD[0]){ fEfficiencyesdTOFPIDD[0]->SetLineColor(3); fEfficiencyesdTOFPIDD[0]->Draw("same"); } TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44); legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency",""); legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p"); legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p"); legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p"); legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p"); legtofeff->Draw("same"); TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500); cefficiencyParamIP->cd(); gStyle->SetOptStat(0); // IP cut efficiencies for(int icentr=0; icentrCalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined"); source = 1; //beauty enhenced if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD); efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined"); if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen); else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD); efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd"); source = 0; //charm mb if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD); efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. charmCombined->Add(mccontainermcD); AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined); efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); source = 1; //beauty mb if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD); efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. beautyCombined->Add(mccontainermcD); AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined); efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD); efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. beautyCombinedesd->Add(mccontaineresdD); AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd); efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); source = 2; //conversion mb if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD); efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. source = 3; //non HFE except for the conversion mb if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD); efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. if(fIPEffCombinedSamples){ fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb } else{ fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only } fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only } if(fBeamType==0){ AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0); fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0); AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0); fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0); } for(int icentr=0; icentrSetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); fipfit->SetLineColor(2); fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R"); fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency"); if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line? else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit"); fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); fipfit->SetLineColor(6); fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R"); fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency"); if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line? else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit"); fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); fipfit->SetLineColor(1); fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R"); fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency"); fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit"); if(fIPParameterizedEff){ fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); fipfitnonhfe->SetLineColor(3); fConversionEff[icentr]->Fit(fipfitnonhfe,"R"); fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency"); fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe"); fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); fipfitnonhfe->SetLineColor(4); fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R"); fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency"); fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe"); } } // draw (for PbPb, only 1st bin) fEfficiencyCharmSigD[0]->SetMarkerStyle(21); fEfficiencyCharmSigD[0]->SetMarkerColor(1); fEfficiencyCharmSigD[0]->SetLineColor(1); fEfficiencyBeautySigD[0]->SetMarkerStyle(21); fEfficiencyBeautySigD[0]->SetMarkerColor(2); fEfficiencyBeautySigD[0]->SetLineColor(2); fEfficiencyBeautySigesdD[0]->SetStats(0); fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21); fEfficiencyBeautySigesdD[0]->SetMarkerColor(6); fEfficiencyBeautySigesdD[0]->SetLineColor(6); fCharmEff[0]->SetMarkerStyle(24); fCharmEff[0]->SetMarkerColor(1); fCharmEff[0]->SetLineColor(1); fBeautyEff[0]->SetMarkerStyle(24); fBeautyEff[0]->SetMarkerColor(2); fBeautyEff[0]->SetLineColor(2); fConversionEff[0]->SetMarkerStyle(24); fConversionEff[0]->SetMarkerColor(3); fConversionEff[0]->SetLineColor(3); fNonHFEEff[0]->SetMarkerStyle(24); fNonHFEEff[0]->SetMarkerColor(4); fNonHFEEff[0]->SetLineColor(4); fEfficiencyCharmSigD[0]->Draw(); fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9); fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5); fEfficiencyBeautySigD[0]->Draw("same"); fEfficiencyBeautySigesdD[0]->Draw("same"); //fCharmEff[0]->Draw("same"); //fBeautyEff[0]->Draw("same"); if(fBeamType==0){ fConversionEffbgc->SetMarkerStyle(25); fConversionEffbgc->SetMarkerColor(3); fConversionEffbgc->SetLineColor(3); fNonHFEEffbgc->SetMarkerStyle(25); fNonHFEEffbgc->SetMarkerColor(4); fNonHFEEffbgc->SetLineColor(4); fConversionEffbgc->Draw("same"); fNonHFEEffbgc->Draw("same"); } else{ fConversionEff[0]->Draw("same"); fNonHFEEff[0]->Draw("same"); } if(fEfficiencyIPBeautyD[0]) fEfficiencyIPBeautyD[0]->Draw("same"); if(fEfficiencyIPBeautyesdD[0]) fEfficiencyIPBeautyesdD[0]->Draw("same"); if( fEfficiencyIPCharmD[0]) fEfficiencyIPCharmD[0]->Draw("same"); if(fIPParameterizedEff){ fEfficiencyIPConversionD[0]->Draw("same"); fEfficiencyIPNonhfeD[0]->Draw("same"); } TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39); legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency",""); legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p"); legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p"); legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p"); legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p"); legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p"); //legipeff->AddEntry(fConversionEff[0],"conversion e","p"); //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p"); legipeff->Draw("same"); gPad->SetGrid(); //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps"); } //____________________________________________________________________________ THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){ // // Return beauty electron IP cut efficiency // const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning const Int_t nCentralitybinning=11;//number of centrality bins 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 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; //Int_t ptpr = 0; Int_t nDim=1; //dimensions of the efficiency weighting grid if(fBeamType==1) { //ptpr=1; nDim=2; //dimensions of the efficiency weighting grid } Int_t nBin[1] = {nPtbinning1}; Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1}; THnSparseF *ipcut; if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin); else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb); for(Int_t idim = 0; idim < nDim; idim++) { if(nDim==1) ipcut->SetBinEdges(idim, kPtRange); if(nDim==2) { ipcut->SetBinEdges(0, kCentralityRange); ipcut->SetBinEdges(1, kPtRange); } } Double_t pt[1]; Double_t weight[nCentralitybinning]; Double_t weightErr[nCentralitybinning]; Double_t contents[2]; for(Int_t a=0;a<11;a++) { weight[a] = 1.0; weightErr[a] = 1.0; } Int_t looppt=nBin[0]; Int_t loopcentr=1; Int_t ibin[2]; if(fBeamType==1) { loopcentr=nBinPbPb[0]; } for(int icentr=0; icentrEval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr); weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1); } } else{ if(fEfficiencyIPBeautyesdD[icentr]){ weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr); weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1); weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1); } } if(fBeamType==1){ contents[0]=icentr; contents[1]=pt[0]; ibin[0]=icentr; ibin[1]=i+1; } if(fBeamType==0){ contents[0]=pt[0]; ibin[0]=i+1; } ipcut->Fill(contents,weight[icentr]); ipcut->SetBinError(ibin,weightErr[icentr]); } } Int_t nDimSparse = ipcut->GetNdimensions(); Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); nBins *= binsvar[iVar]; } Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) // loop that sets 0 error in each bin for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere } delete[] binsvar; delete[] binfill; return ipcut; } //____________________________________________________________________________ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){ // // Return PID x IP cut efficiency // const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning const Int_t nCentralitybinning=11;//number of centrality bins 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 Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; //Int_t ptpr = 0; Int_t nDim=1; //dimensions of the efficiency weighting grid if(fBeamType==1) { //ptpr=1; nDim=2; //dimensions of the efficiency weighting grid } Int_t nBin[1] = {nPtbinning1}; Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1}; THnSparseF *pideff; if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin); else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb); for(Int_t idim = 0; idim < nDim; idim++) { if(nDim==1) pideff->SetBinEdges(idim, kPtRange); if(nDim==2) { pideff->SetBinEdges(0, kCentralityRange); pideff->SetBinEdges(1, kPtRange); } } Double_t pt[1]; Double_t weight[nCentralitybinning]; Double_t weightErr[nCentralitybinning]; Double_t contents[2]; for(Int_t a=0;aEval(0); // assume we have constant TRD+TPC PID efficiency for(int i=0; iEval(pt[0]); else{ printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr); } if(fEfficiencyesdTOFPIDD[icentr]) tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]); else{ printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr); } //tof pid eff x tpc pid eff x ip cut eff if(fIPParameterizedEff){ if(source==0) { if(fEfficiencyIPCharmD[icentr]){ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr); weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1); } } else if(source==2) { if(fEfficiencyIPConversionD[icentr]){ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr); weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1); } } else if(source==3) { if(fEfficiencyIPNonhfeD[icentr]){ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr); weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1); } } } else{ if(source==0){ if(fEfficiencyIPCharmD[icentr]){ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]); weightErr[icentr] = 0; } else{ printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr); weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1); } } else if(source==2){ if(fBeamType==0){ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1); } else{ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1); } } else if(source==3){ if(fBeamType==0){ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1); } else{ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1); } } } if(fBeamType==1){ contents[0]=icentr; contents[1]=pt[0]; ibin[0]=icentr; ibin[1]=i+1; } if(fBeamType==0){ contents[0]=pt[0]; ibin[0]=i+1; } pideff->Fill(contents,weight[icentr]); pideff->SetBinError(ibin,weightErr[icentr]); } } Int_t nDimSparse = pideff->GetNdimensions(); Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); nBins *= binsvar[iVar]; } Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) // loop that sets 0 error in each bin for (Long_t iBin=0; iBinSetBinEdges(0, kPtRange); if(fBeamType==1) { // brawspectra->SetBinEdges(0, centralityBins); brawspectra->SetBinEdges(0, kCentralityRange); brawspectra->SetBinEdges(1, kPtRange); } Double_t pt[1]; Double_t yields= 0.; Double_t valuesb[2]; //Int_t looppt=nBin[0]; Int_t loopcentr=1; if(fBeamType==1) { loopcentr=nBinPbPb[0]; } for(int icentr=0; icentrGetNbinsX(); i++){ pt[0]=(kPtRange[i]+kPtRange[i+1])/2.; yields = fBSpectrum2ndMethod->GetBinContent(i+1); if(fBeamType==1) { valuesb[0]=icentr; valuesb[1]=pt[0]; } if(fBeamType==0) valuesb[0]=pt[0]; brawspectra->Fill(valuesb,yields); } } Int_t nDimSparse = brawspectra->GetNdimensions(); Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); nBins *= binsvar[iVar]; } Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) // loop that sets 0 error in each bin for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere } rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr); new TCanvas; fBSpectrum2ndMethod->SetMarkerStyle(24); fBSpectrum2ndMethod->Draw("p"); hRawBeautySpectra->SetMarkerStyle(25); hRawBeautySpectra->Draw("samep"); delete[] binfill; delete[] binsvar; return rawBeautyContainer; } //__________________________________________________________________________ void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){ // // Calculate non HFE sys // // if(!fNonHFEsyst) return; Double_t evtnorm[1] = {0.0}; if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels]; AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels]; AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors AliCFDataGrid *bgNonHFEGrid[kBgLevels]; AliCFDataGrid *bgConvGrid[kBgLevels]; Int_t stepbackground = 3; Int_t* bins=new Int_t[1]; const Char_t *bgBase[2] = {"pi0","eta"}; bins[0]=fConversionEff[centrality]->GetNbinsX(); AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins); AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins); for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){ 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); weightedConversionContainer->SetGrid(GetPIDxIPEff(2)); convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0); 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); weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3)); nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0); } bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone(); for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){ bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]); } if(!fEtaSyst) bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]); bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone(); 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) bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]); } if(!fEtaSyst) bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]); bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone(); bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]); if(fEtaSyst){ bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]); } } //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) AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone(); bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.); bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone(); bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.); //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0); hBaseErrors[iErr][0]->Scale(-1.); hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr])); hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0); if(!fEtaSyst)break; } //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 TH1D *hSpeciesErrors[kElecBgSources-1]; for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ if(fEtaSyst && (iSource == 1))continue; hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0); TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0); hSpeciesErrors[iSource-1]->Add(hNonHFEtemp); hSpeciesErrors[iSource-1]->Scale(0.3); } //Int_t firstBgSource = 0;//if eta systematics are not from scaling //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined! TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone(); TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone(); TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone(); TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone(); TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone(); for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){ Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp; pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin); pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin); if(fEtaSyst){ etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin); etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin); } else{ etaErrLow = etaErrUp = 0.;} Double_t sqrsumErrs= 0; for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ if(fEtaSyst && (iSource == 1))continue; Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin); sqrsumErrs+=(scalingErr*scalingErr); } for(Int_t iErr = 0; iErr < 2; iErr++){ for(Int_t iLevel = 0; iLevel < 2; iLevel++){ hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin)); } if(!fEtaSyst)break; } hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs)); hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs)); hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin)); hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin)); hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin)); } TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600); cPiErrors->cd(); cPiErrors->SetLogx(); cPiErrors->SetLogy(); hBaseErrors[0][0]->Draw(); //hBaseErrors[0][1]->SetMarkerColor(kBlack); //hBaseErrors[0][1]->SetLineColor(kBlack); //hBaseErrors[0][1]->Draw("SAME"); if(fEtaSyst){ hBaseErrors[1][0]->Draw("SAME"); hBaseErrors[1][0]->SetMarkerColor(kBlack); hBaseErrors[1][0]->SetLineColor(kBlack); //hBaseErrors[1][1]->SetMarkerColor(13); //hBaseErrors[1][1]->SetLineColor(13); //hBaseErrors[1][1]->Draw("SAME"); } //hOverallBinScaledErrsUp->SetMarkerColor(kBlue); //hOverallBinScaledErrsUp->SetLineColor(kBlue); //hOverallBinScaledErrsUp->Draw("SAME"); hOverallBinScaledErrsLow->SetMarkerColor(kGreen); hOverallBinScaledErrsLow->SetLineColor(kGreen); hOverallBinScaledErrsLow->Draw("SAME"); hScalingErrors->SetLineColor(kBlue); hScalingErrors->Draw("SAME"); TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95); lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error"); //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error"); if(fEtaSyst){ lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error"); //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error"); } lPiErr->AddEntry(hScalingErrors, "scaling error"); lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics"); //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics"); lPiErr->Draw("SAME"); //Normalize errors TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone(); TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone(); hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!) hLowSystScaled->Scale(evtnorm[0]); TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone(); TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone(); //histograms to be normalized to TGraphErrors CorrectFromTheWidth(hNormAllSystErrUp); CorrectFromTheWidth(hNormAllSystErrLow); TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs"); cNormOvErrs->cd(); cNormOvErrs->SetLogx(); cNormOvErrs->SetLogy(); TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp); TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow); gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors"); gOverallSystErrUp->SetMarkerColor(kBlack); gOverallSystErrUp->SetLineColor(kBlack); gOverallSystErrLow->SetMarkerColor(kRed); gOverallSystErrLow->SetLineColor(kRed); gOverallSystErrUp->Draw("AP"); gOverallSystErrLow->Draw("PSAME"); TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89); lAllSys->AddEntry(gOverallSystErrLow,"lower","p"); lAllSys->AddEntry(gOverallSystErrUp,"upper","p"); lAllSys->Draw("same"); AliCFDataGrid *bgYieldGrid; if(fEtaSyst){ 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. } bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone(); TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0); TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone(); hRelErrUp->Divide(hBgYield); TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone(); hRelErrLow->Divide(hBgYield); TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs"); cRelErrs->cd(); cRelErrs->SetLogx(); hRelErrUp->SetTitle("Relative error of non-HFE background yield"); hRelErrUp->Draw(); hRelErrLow->SetLineColor(kBlack); hRelErrLow->Draw("SAME"); TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95); lRel->AddEntry(hRelErrUp, "upper"); lRel->AddEntry(hRelErrLow, "lower"); lRel->Draw("SAME"); //CorrectFromTheWidth(hBgYield); //hBgYield->Scale(evtnorm[0]); //write histograms/TGraphs to file TFile *output = new TFile("systHists.root","recreate"); hBgYield->SetName("hBgYield"); hBgYield->Write(); hRelErrUp->SetName("hRelErrUp"); hRelErrUp->Write(); hRelErrLow->SetName("hRelErrLow"); hRelErrLow->Write(); hUpSystScaled->SetName("hOverallSystErrUp"); hUpSystScaled->Write(); hLowSystScaled->SetName("hOverallSystErrLow"); hLowSystScaled->Write(); gOverallSystErrUp->SetName("gOverallSystErrUp"); gOverallSystErrUp->Write(); gOverallSystErrLow->SetName("gOverallSystErrLow"); gOverallSystErrLow->Write(); output->Close(); delete output; } //____________________________________________________________ void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){ // // Unfold backgrounds to check its sanity // AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC); //AliCFContainer *mcContainer = GetContainer(kMCContainerMC); if(!mcContainer){ AliError("MC Container not available"); return; } if(!fCorrelation){ AliError("No Correlation map available"); return; } // Data AliCFDataGrid *dataGrid = 0x0; dataGrid = bgsubpectrum; // Guessed AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding); THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid(); // Efficiency AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue); // Unfold background spectra Int_t nDim=1; if(fBeamType==0)nDim = 1; if(fBeamType==1)nDim = 2; AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse); if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors(); unfolding.SetMaxNumberOfIterations(fNumberOfIterations); if(fSetSmoothing) unfolding.UseSmoothing(); unfolding.Unfold(); // Results THnSparse* result = unfolding.GetUnfolded(); TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600); if(fBeamType==1) { ctest->Divide(2,2); ctest->cd(1); result->GetAxis(0)->SetRange(1,1); TH1D* htest1=(TH1D*)result->Projection(0); htest1->Draw(); ctest->cd(2); result->GetAxis(0)->SetRange(1,1); TH1D* htest2=(TH1D*)result->Projection(1); htest2->Draw(); ctest->cd(3); result->GetAxis(0)->SetRange(6,6); TH1D* htest3=(TH1D*)result->Projection(0); htest3->Draw(); ctest->cd(4); result->GetAxis(0)->SetRange(6,6); TH1D* htest4=(TH1D*)result->Projection(1); htest4->Draw(); } TGraphErrors* unfoldedbgspectrumD = Normalize(result); if(!unfoldedbgspectrumD) { AliError("Unfolded background spectrum doesn't exist"); } else{ TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate"); if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum"); if(fBeamType==1) { Int_t centr=1; result->GetAxis(0)->SetRange(centr,centr); unfoldedbgspectrumD = Normalize(result,centr-1); unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20"); centr=6; result->GetAxis(0)->SetRange(centr,centr); unfoldedbgspectrumD = Normalize(result,centr-1); unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80"); } file->Close(); } }