fAnalyzedEvents(0) ,
fCentralityId(0) ,
fCentStart(0.),
- fCentStop(0.)
+ fCentStop(0.),
+ fHistHBTbefore(NULL),
+ fHistHBTafter(NULL),
+ fHistConversionbefore(NULL),
+ fHistConversionafter(NULL)
{
// Default constructor
fHistNN[i] = NULL;
}
-
- //QA histograms
- fHistHBTbefore = NULL;
- fHistHBTafter = NULL;
- fHistConversionbefore = NULL;
- fHistConversionafter = NULL;
-
}
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 += kBFAnalysisType[iAnalysisType];
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::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: "<<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;
+ 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,Bool_t correctWithEfficiency) {
+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){
- fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/effFromConvolutionAllCent.root");
+ 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;
// - 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){
+ if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency && !correctWithMixed){
- 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* 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;
// - 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){
+ if(iAnalysisType == kPhi && correctWithEfficiency && !correctWithMixed){
- TH1F* hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
- TH1F* hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
+ 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;
}
}
+ // 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)) {
}
// do the acceptance correction (only for Eta and etaWindow > 0)
- if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency){
+ if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency && !correctWithMixed){
for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff