X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGCF%2FEBYE%2FBalanceFunctions%2FAliBalance.cxx;h=6c8d7414863ab52054c70f77d183bfd4e7651f5e;hb=da8b2d30c9fba2b21d40e79c7f82b62f0526f94b;hp=f0b4f9a8144714e4c617030db382909fbf797418;hpb=e1d6ee8e877d981b72347486668704ac7b20cf76;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalance.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalance.cxx index f0b4f9a8144..6c8d7414863 100644 --- a/PWGCF/EBYE/BalanceFunctions/AliBalance.cxx +++ b/PWGCF/EBYE/BalanceFunctions/AliBalance.cxx @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include #include @@ -37,19 +39,27 @@ #include "AliBalance.h" +using std::cout; +using std::cerr; +using std::endl; + ClassImp(AliBalance) //____________________________________________________________________// AliBalance::AliBalance() : TObject(), - bShuffle(kFALSE), - bHBTcut(kFALSE), - bConversionCut(kFALSE), + fShuffle(kFALSE), + fHBTcut(kFALSE), + fConversionCut(kFALSE), fAnalysisLevel("ESD"), fAnalyzedEvents(0) , fCentralityId(0) , fCentStart(0.), - fCentStop(0.) + fCentStop(0.), + fHistHBTbefore(NULL), + fHistHBTafter(NULL), + fHistConversionbefore(NULL), + fHistConversionafter(NULL) { // Default constructor @@ -99,14 +109,18 @@ AliBalance::AliBalance() : //____________________________________________________________________// AliBalance::AliBalance(const AliBalance& balance): TObject(balance), - bShuffle(balance.bShuffle), - bHBTcut(balance.bHBTcut), - bConversionCut(balance.bConversionCut), + fShuffle(balance.fShuffle), + fHBTcut(balance.fHBTcut), + fConversionCut(balance.fConversionCut), fAnalysisLevel(balance.fAnalysisLevel), fAnalyzedEvents(balance.fAnalyzedEvents), fCentralityId(balance.fCentralityId), fCentStart(balance.fCentStart), - fCentStop(balance.fCentStop) { + fCentStop(balance.fCentStop), + fHistHBTbefore(balance.fHistHBTbefore), + fHistHBTafter(balance.fHistHBTafter), + fHistConversionbefore(balance.fHistConversionbefore), + fHistConversionafter(balance.fHistConversionafter) { //copy constructor for(Int_t i = 0; i < ANALYSIS_TYPES; i++){ fNn[i] = balance.fNn[i]; @@ -190,42 +204,58 @@ void AliBalance::SetInterval(Int_t iAnalysisType, //____________________________________________________________________// void AliBalance::InitHistograms() { //Initialize the histograms + + // global switch disabling the reference + // (to avoid "Replacing existing TH1" if several wagons are created in train) + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + TString histName; for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) { - histName = "fHistP"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistP"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]); - histName = "fHistN"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistN"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]); - histName = "fHistPN"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistPN"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]); - histName = "fHistNP"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistNP"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]); - histName = "fHistPP"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistPP"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]); - histName = "fHistNN"; histName += gBFAnalysisType[iAnalysisType]; - if(bShuffle) histName.Append("_shuffle"); + histName = "fHistNN"; histName += kBFAnalysisType[iAnalysisType]; + if(fShuffle) histName.Append("_shuffle"); if(fCentralityId) histName += fCentralityId.Data(); fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]); } + + // QA histograms + fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200); + fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200); + fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200); + fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200); + + TH1::AddDirectory(oldStatus); + } //____________________________________________________________________// void AliBalance::PrintAnalysisSettings() { + //prints the analysis settings Printf("======================================"); Printf("Analysis level: %s",fAnalysisLevel.Data()); @@ -394,13 +424,15 @@ void AliBalance::CalculateBalance(Float_t fCentrality,vector **chargeV if(dphi>180) dphi = 360 - dphi; //dphi should be between 0 and 180! // HBT like cut - if(bHBTcut && charge1 * charge2 > 0){ + if(fHBTcut && charge1 * charge2 > 0){ //if( dphi < 3 || deta < 0.01 ){ // VERSION 1 // continue; // VERSION 2 (Taken from DPhiCorrelations) // the variables & cuthave been developed by the HBT group // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700 + + fHistHBTbefore->Fill(deta,dphi); // optimization if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations] @@ -440,12 +472,16 @@ void AliBalance::CalculateBalance(Float_t fCentrality,vector **chargeV } } } + fHistHBTafter->Fill(deta,dphi); } // conversions - if(bConversionCut){ + if(fConversionCut){ if (charge1 * charge2 < 0) { + + fHistConversionbefore->Fill(deta,dphi); + Float_t m0 = 0.510e-3; Float_t tantheta1 = 1e10; @@ -469,6 +505,7 @@ void AliBalance::CalculateBalance(Float_t fCentrality,vector **chargeV //AliInfo(Form("Conversion: Removed track pair %d %d with [[%f %f] %f %f] %d %d <- %f %f %f %f %f %f ", i, j, deta, dphi, masssqu, charge1, charge2,eta1,eta2,phi1,phi2,pt1,pt2)); continue; } + fHistConversionafter->Fill(deta,dphi); } } @@ -746,19 +783,18 @@ TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) { //____________________________________________________________________// void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) { //Prints the calculated width of the BF and its error - Double_t x[MAXIMUM_NUMBER_OF_STEPS]; Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0; Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0; Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0; Double_t deltaBalP2 = 0.0, integral = 0.0; Double_t deltaErrorNew = 0.0; - cout<<"=================================================="<GetBinContent(i)<<"\t Error: "<GetBinError(i)<<"\t bin: "<GetBinCenter(i)<GetBinContent(i)<<"\t Error: "<GetBinError(i)<<"\t bin: "<GetBinCenter(i)<GetBinCenter(i); gSumBi += gHistBalance->GetBinContent(i); @@ -775,18 +811,20 @@ void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) { deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2); Double_t integralError = TMath::Sqrt(deltaBalP2); - - Double_t delta = gSumBiXi / gSumBi; + integralError *= 1.0; + + Double_t delta = gSumBiXi / gSumBi; delta *= 1.0; Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) ); - cout<<"Analysis type: "<(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX)); TH1D *hTemp6 = dynamic_cast(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX)); + // get the file with the efficiency matrices + // withAcceptanceOnly: Data single distributions are normalized to 1 (efficiency not taken into account) + // else : Data single distributions are normalized to give single particle efficiency of MC + TFile *fEfficiencyMatrix = NULL; + if(correctWithEfficiency || correctWithMixed){ + if(correctWithAcceptanceOnly) fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/accOnlyFromConvolutionAllCent.root"); + else fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/effFromConvolutionAllCent.root"); + if(!fEfficiencyMatrix){ + AliError("Efficiency histogram file not found"); + return NULL; + } + } + + // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations) + // - single particle efficiencies from MC (AliAnalysiTaskEfficiency) + // - two particle efficiencies from convolution of data single particle distributions + // (normalized to single particle efficiency) + if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency && !correctWithMixed){ + + TH1F* hEffP = NULL; + TH1F* hEffN = NULL; + TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + + // take the data distributions + if(correctWithAcceptanceOnly){ + hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + } + // take the MC distributions + else{ + hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax)); + hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax)); + } + + if( !hEffP || !hEffN || !hEffPP || !hEffNN || !hEffPN){ + AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + return NULL; + } + + for(Int_t iBin = 0; iBin < hEffP->GetNbinsX(); iBin++){ + hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1)))); + hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1)))); + + hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1)))); + hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1)))); + } + + for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){ + + hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1)))); + hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1)))); + hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1)))); + hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1)))); + hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1)))); + hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1)))); + hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1)))); + hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1)))); + + } + + // TF1 *fPP = new TF1("fPP","pol1",0,1.6); // phase space factor + efficiency for ++ + // fPP->SetParameters(0.736466,-0.461529); + // TF1 *fNN = new TF1("fNN","pol1",0,1.6); // phase space factor + efficiency for -- + // fNN->SetParameters(0.718616,-0.450473); + // TF1 *fPN = new TF1("fPN","pol1",0,1.6); // phase space factor + efficiency for +- + // fPN->SetParameters(0.727507,-0.455981); + + // for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){ + // hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1))); + // hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1))); + // } + } + + // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations) + // - single particle efficiencies from MC (AliAnalysiTaskEfficiency) + // - two particle efficiencies from convolution of data single particle distributions + // (normalized to single particle efficiency) + if(iAnalysisType == kPhi && correctWithEfficiency && !correctWithMixed){ + + TH1F* hEffPhiP = NULL; + TH1F* hEffPhiN = NULL; + TH1F* hEffPhiPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffPhiNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffPhiPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + + // take the data distributions + if(correctWithAcceptanceOnly){ + hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + } + // take the MC distributions + else{ + hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_MC",centrMin,centrMax)); + hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_MC",centrMin,centrMax)); + } + + if( !hEffPhiP || !hEffPhiN || !hEffPhiPP || !hEffPhiNN || !hEffPhiPN){ + AliError("Efficiency (phi) histograms not found"); + return NULL; + } + + for(Int_t iBin = 0; iBin < hEffPhiP->GetNbinsX(); iBin++){ + hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1)))); + hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1)))); + + hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1)))); + hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1)))); + } + + for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){ + + hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1)))); + hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1)))); + hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1)))); + hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1)))); + hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1)))); + hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1)))); + hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1)))); + hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1)))); + + } + } + + // do the correction with the event mixing directly! + if(correctWithMixed){ + + // take the MC distributions (for average efficiency) + TH1F* hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax)); + TH1F* hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax)); + + TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax)); + + if( !hEffP || !hEffN){ + AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax)); + return NULL; + } + + if(hMixed[0] && hMixed[1] && hMixed[2] && hMixed[3]){ + + // scale to average efficiency in the pt region (0.3-1.5) and |eta| < 0.8 + // by multiplying the average single particle efficiencies from HIJING + // here we assume that the distributions are 1: + // - in the integral for dphi (for averaging over sector structure) + // - in the maximum for deta + Double_t normPMC = (Double_t)hEffP->Integral()/(Double_t)hEffP->GetNbinsX(); + Double_t normNMC = (Double_t)hEffN->Integral()/(Double_t)hEffN->GetNbinsX(); + Double_t normPPMC = (Double_t)hEffPP->Integral()/(Double_t)hEffPP->GetNbinsX(); + Double_t normNNMC = (Double_t)hEffNN->Integral()/(Double_t)hEffNN->GetNbinsX(); + Double_t normPNMC = (Double_t)hEffPN->Integral()/(Double_t)hEffPN->GetNbinsX(); + + hMixed[0]->Scale(normPNMC); + hMixed[1]->Scale(normPNMC); + hMixed[2]->Scale(normNNMC); + hMixed[3]->Scale(normPPMC); + + // divide by event mixing + hTemp1->Divide(hMixed[0]); + hTemp2->Divide(hMixed[1]); + hTemp3->Divide(hMixed[2]); + hTemp4->Divide(hMixed[3]); + + // scale also single histograms with average efficiency + hTemp5->Scale(1./normNMC); + hTemp6->Scale(1./normPMC); + + } + else{ + AliError("Correction with EventMixing requested, but not all Histograms there!"); + return NULL; + } + } + + if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) { hTemp1->Sumw2(); hTemp2->Sumw2(); hTemp3->Sumw2(); hTemp4->Sumw2(); hTemp1->Add(hTemp3,-2.); - hTemp1->Scale(1./hTemp5->GetEntries()); + hTemp1->Scale(1./hTemp5->Integral()); hTemp2->Add(hTemp4,-2.); - hTemp2->Scale(1./hTemp6->GetEntries()); + hTemp2->Scale(1./hTemp6->Integral()); gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.); gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]); } // do the acceptance correction (only for Eta and etaWindow > 0) - if(iAnalysisType == kEta && etaWindow > 0){ + if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency && !correctWithMixed){ for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){ Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1); Double_t corrected = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow ); gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected); + gHistBalanceFunctionHistogram->SetBinError(iBin+1,corrected/notCorrected*gHistBalanceFunctionHistogram->GetBinError(iBin+1)); } } + if(fEfficiencyMatrix) fEfficiencyMatrix->Close(); + PrintResults(iAnalysisType,gHistBalanceFunctionHistogram); return gHistBalanceFunctionHistogram;