#include <Riostream.h>
#include <TMath.h>
#include <TAxis.h>
+#include <TFile.h>
+#include <TF1.h>
#include <TH2D.h>
#include <TLorentzVector.h>
#include <TObjArray.h>
#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
//____________________________________________________________________//
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];
//____________________________________________________________________//
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());
px2 = chargeVector[4]->at(j);
py2 = chargeVector[5]->at(j);
pz2 = chargeVector[6]->at(j);
- pt2 = chargeVector[7]->at(i);
+ pt2 = chargeVector[7]->at(j);
energy2 = chargeVector[8]->at(j);
-
+
// filling the arrays
// RAPIDITY
if(dphi>180) dphi = 360 - dphi; //dphi should be between 0 and 180!
// HBT like cut
- if(bHBTcut){
+ 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]
{
+
+ // phi in rad
+ Float_t phi1rad = phi1*TMath::DegToRad();
+ Float_t phi2rad = phi2*TMath::DegToRad();
+
// check first boundaries to see if is worth to loop and find the minimum
- Float_t dphistar1 = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, 0.8, bSign);
- Float_t dphistar2 = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, 2.5, bSign);
+ Float_t dphistar1 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 0.8, bSign);
+ Float_t dphistar2 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 2.5, bSign);
const Float_t kLimit = 0.02 * 3;
Float_t dphistarminabs = 1e5;
Float_t dphistarmin = 1e5;
- if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
+
+ if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
{
for (Double_t rad=0.8; rad<2.51; rad+=0.01)
{
- Float_t dphistar = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, rad, bSign);
-
+ Float_t dphistar = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, rad, bSign);
Float_t dphistarabs = TMath::Abs(dphistar);
if (dphistarabs < dphistarminabs)
dphistarminabs = dphistarabs;
}
}
-
+
if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
{
- //Printf("Removed track pair %d %d with %f %f %f %f %d %f %f %d %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
+ //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
continue;
}
}
}
+ 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;
+
+ // phi in rad
+ Float_t phi1rad = phi1*TMath::DegToRad();
+ Float_t phi2rad = phi2*TMath::DegToRad();
if (eta1 < -1e-10 || eta1 > 1e-10)
tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
-
- Float_t mass = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
-
- if (mass < 0.04*0.04){
- //Printf("Removed track pair %d %d with %f %f %d %d ", i, j, deta, mass, charge1, charge2);
+
+ Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
+
+ if (masssqu < 0.04*0.04){
+ //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);
}
}
//____________________________________________________________________//
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<<"=================================================="<<endl;
- for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
- x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
- //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
- }
- //cout<<"=================================================="<<endl;
+ // cout<<"=================================================="<<endl;
+ // for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
+ // x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
+ // cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+ // }
+ // cout<<"=================================================="<<endl;
for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
gSumXi += gHistBalance->GetBinCenter(i);
gSumBi += gHistBalance->GetBinContent(i);
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: "<<gBFAnalysisType[iAnalysisType].Data()<<endl;
- cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
- cout<<"New error: "<<deltaErrorNew<<endl;
- cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
- cout<<"=================================================="<<endl;
+ deltaError *= 1.0;
+ // cout<<"Analysis type: "<<kBFAnalysisType[iAnalysisType].Data()<<endl;
+ // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+ // cout<<"New error: "<<deltaErrorNew<<endl;
+ // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+ // cout<<"=================================================="<<endl;
}
//____________________________________________________________________//
-TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow) {
+TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly, Bool_t correctWithMixed, TH1D *hMixed[4]) {
//Returns the BF histogram, extracted from the 6 TH2D objects
//(private members) of the AliBalance class.
//
TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
TH1D *hTemp6 = dynamic_cast<TH1D *>(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;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff