+//__________________________________________________________________________
+void AliHFEspectrum::CalculateNonHFEsyst(){
+ //
+ //Calculate systematic errors for non-HFE and conversion background,
+ //based on correlated errors from pi0 input and uncorrelated errors
+ //from m_t scaling
+ //
+ if(!fNonHFEsyst)
+ return;
+
+ Double_t evtnorm[1] = {0.0};
+ if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+
+ AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
+ AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
+
+ AliCFDataGrid *bgLevelGrid[kBgLevels];
+ AliCFDataGrid *bgNonHFEGrid[kBgLevels];
+ AliCFDataGrid *bgConvGrid[kBgLevels];
+
+ Int_t stepbackground = 3;
+ Int_t* bins=new Int_t[1];
+ bins[0]=fConversionEff->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",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
+ weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
+ convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
+
+ nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
+ weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
+ nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
+ }
+
+ bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
+ }
+
+ bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays
+ bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
+ }
+
+ bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
+ bgLevelGrid[iLevel]->Add(bgNonHFEGrid[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)
+ AliCFDataGrid *bgErrorGrid[2];
+ bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone();
+ bgErrorGrid[1] = (AliCFDataGrid*)bgLevelGrid[2]->Clone();
+ bgErrorGrid[0]->Add(bgLevelGrid[0],-1.);
+ bgErrorGrid[1]->Add(bgLevelGrid[0],-1.);
+
+ //plot absolute differences between limit yields (upper/lower limit, based on pi0 errors) and best estimate
+ TH1D* hpiErrors[2];
+ hpiErrors[0] = (TH1D*)bgErrorGrid[0]->Project(0);
+ hpiErrors[0]->Scale(-1.);
+ hpiErrors[0]->SetTitle("Absolute systematic errors from non-HF meson decays and conversions");
+ hpiErrors[1] = (TH1D*)bgErrorGrid[1]->Project(0);
+
+
+
+ //Calculate the scaling errors for electrons from all mesons except for pions: 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++){
+ 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);
+ }
+
+ TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[0]->Clone();
+ TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[0]->Clone();
+ TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[0]->Clone();
+
+ TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
+ TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
+
+ for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
+ Double_t pi0basedErrLow = hpiErrors[0]->GetBinContent(iBin);
+ Double_t pi0basedErrUp = hpiErrors[1]->GetBinContent(iBin);
+
+ Double_t sqrsumErrs= 0;
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
+ sqrsumErrs+=(scalingErr*scalingErr);
+ }
+ for(Int_t iErr = 0; iErr < 2; iErr++){
+ hpiErrors[iErr]->SetBinContent(iBin,hpiErrors[iErr]->GetBinContent(iBin)/hpiErrors[iErr]->GetBinWidth(iBin));
+ }
+ hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+sqrsumErrs));
+ hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+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));
+ }
+
+
+ // /hOverallSystErrUp->GetBinWidth(iBin))
+
+ TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
+ cPiErrors->cd();
+ cPiErrors->SetLogx();
+ cPiErrors->SetLogy();
+ hpiErrors[0]->Draw();
+ hpiErrors[1]->SetMarkerColor(kBlack);
+ hpiErrors[1]->SetLineColor(kBlack);
+ hpiErrors[1]->Draw("SAME");
+ hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
+ hOverallBinScaledErrsUp->SetLineColor(kBlue);
+ hOverallBinScaledErrsUp->Draw("SAME");
+ hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
+ hOverallBinScaledErrsLow->SetLineColor(kGreen);
+ hOverallBinScaledErrsLow->Draw("SAME");
+ hScalingErrors->SetLineColor(11);
+ hScalingErrors->Draw("SAME");
+
+ TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
+ lPiErr->AddEntry(hpiErrors[0],"Lower error from pion error");
+ lPiErr->AddEntry(hpiErrors[1],"Upper error from pion 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;
+ bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[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;
+}