From 7827f01edd54f44301e4b51f526102065e379ac9 Mon Sep 17 00:00:00 2001 From: miweber Date: Mon, 30 Jun 2014 20:42:55 +0200 Subject: [PATCH] Adding FWHM to determine BF width (Alis Rodriguez Manso ) --- PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx | 223 +++++++++++++++++- PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h | 5 + PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C | 77 ++++-- 3 files changed, 281 insertions(+), 24 deletions(-) diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx index 011c640fbeb..f93257604e3 100644 --- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx +++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include "AliVParticle.h" #include "AliMCParticle.h" @@ -2431,7 +2433,6 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin, } //____________________________________________________________________// - Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM, Double_t &mean, Double_t &meanError, Double_t &sigma, Double_t &sigmaError, @@ -2654,3 +2655,223 @@ Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, AliFatal(Form("Tag %s not found in %s", tag, configuration)); return 0; } + +//____________________________________________________________________// +Double_t AliBalancePsi::GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist, + Double_t &fwhm_spline, Double_t &fwhmError) { + + // helper method to calculate the fwhm and errors of a TH1D + if(gHist){ + Int_t repeat = 10000; + + if (gDeltaEtaPhi == 1){ + fwhm_spline = 0.; + fwhmError = 0.; + gHist->Smooth(4); + + Double_t xmin = gHist->GetXaxis()->GetXmin(); + Double_t xmax = gHist->GetXaxis()->GetXmax(); + + //+++++++FWHM TSpline+++++++++++++++++++++++++++// + TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); + spline3->SetLineColor(kGreen+2); + spline3->SetLineWidth(2); + //spline3->Draw("SAME"); + + Int_t Nbin = gHist->GetNbinsX(); + Int_t bin_max_hist_y = gHist->GetMaximumBin(); + Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y); + Double_t max_spline = spline3->Eval(bins_center_x); + Double_t y_spline = -999; + Double_t y_spline_r = -999; + Double_t x_spline_l = -999; + Double_t x_spline_r = -999; + + for (Int_t i= 0; i < bin_max_hist_y*1000; i++){ + y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000)); + if (y_spline > max_spline/2){ + x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000); + break; + } + } + for (Int_t j= 0; j < bin_max_hist_y*1000; j++){ + y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000)); + if (y_spline_r < max_spline/2){ + x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000); + break; + } + } + fwhm_spline = x_spline_r - x_spline_l; + //+++++++FWHM TSpline++++++++++++++++++++++++// + + //+++++++Error Calculation SPLINE++++++++++++// + Double_t dmeans,dsigmas; + TSpline3 *spline1; + Int_t bin_max_hist_y1; + Double_t bins_center_x1; + Double_t max_spline1; + Double_t y_spline1 = - 999.; + Double_t y_spline_r1 = - 999.; + Double_t x_spline_l1 = -999.; + Double_t x_spline_r1 = -999.; + Double_t fwhm_spline1 = 0.; + Double_t fwhm_T = 0.; + + TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax); + TRandom3 *rgauss = new TRandom3(0); + TH1F *hists = new TH1F("hists","hists",1000,1.,3); + + for (Int_t f=0; fGetBinContent(h); + dsigmas = gHist->GetBinError(h); + Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas); + hEmpty2->SetBinContent(h,rgauss_point); + + //++++++++++++// + hEmpty2->Smooth(4); + //++++++++++++// + } + spline1 = new TSpline3(hEmpty2,"b2e2",0,0); + spline1->SetLineColor(kMagenta+1); + spline1->SetLineWidth(1); + //spline1->Draw("SAME"); + + bin_max_hist_y1 = hEmpty2->GetMaximumBin(); + bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1); + max_spline1 = spline1->Eval(bins_center_x1); + + for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){ + y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000)); + if (y_spline1 > max_spline1/2){ + x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000); + break; + } + } + for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){ + y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000)); + if (y_spline_r1 < max_spline1/2){ + x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000); + break; + } + } + + fwhm_spline1 = x_spline_r1 - x_spline_l1; + hists->Fill(fwhm_spline1); + fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1); + } + + fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1))); + //+++++++Error Calculation SPLINE+++++++++++// + } + else{ + fwhm_spline = 0.; + fwhmError = 0.; + gHist->Smooth(4); + + Double_t xmin = gHist->GetXaxis()->GetXmin(); + Double_t xmax = gHist->GetXaxis()->GetXmax(); + + //+++++++FWHM TSpline+++++++++++++++++++++++++++// + TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); + spline3->SetLineColor(kGreen+2); + spline3->SetLineWidth(2); + //spline3->Draw("SAME"); + + Int_t Nbin = gHist->GetNbinsX(); + Int_t bin_max_hist_y = gHist->GetMaximumBin(); + Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y); + Double_t max_spline = spline3->Eval(bins_center_x); + + Int_t bin_min_hist_y = gHist->GetMinimumBin(); + Double_t bins_center_xmin = gHist->GetBinCenter(bin_min_hist_y); + Double_t min_spline = spline3->Eval(bins_center_xmin); + + Double_t y_spline = -999.; + Double_t y_spline_r = -999.; + Double_t x_spline_l = -999.; + Double_t x_spline_r = -999.; + + for (Int_t i= 0; i < bin_max_hist_y*1000; i++){ + y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000)); + if (y_spline > (max_spline+min_spline)/2){ + x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000); + break; + } + } + for (Int_t j= 0; j < bin_max_hist_y*1000; j++){ + y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000)); + if (y_spline_r < (max_spline+min_spline)/2){ + x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000); + break; + } + } + fwhm_spline = x_spline_r - x_spline_l; + //+++++++FWHM TSpline++++++++++++++++++++++++// + + //+++++++Error Calculation SPLINE++++++++++++// + Double_t dmeans,dsigmas; + TSpline3 *spline1; + Int_t bin_max_hist_y1; + Double_t bins_center_x1; + Double_t max_spline1; + Double_t y_spline1 = -999.; + Double_t y_spline_r1 = -999.; + Double_t x_spline_l1 = -999.; + Double_t x_spline_r1 = -999.; + Double_t fwhm_spline1 = 0.; + Double_t fwhm_T = 0.; + + TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax); + TRandom3 *rgauss = new TRandom3(0); + Int_t bin_min_hist_y1; + Double_t bins_center_xmin1,min_spline1; + + for (Int_t f=0; fGetBinContent(h); + dsigmas = gHist->GetBinError(h); + Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas); + hEmpty2->SetBinContent(h,rgauss_point); + + //++++++++++++// + hEmpty2->Smooth(4); + //++++++++++++// + } + spline1 = new TSpline3(hEmpty2,"b2e2",0,0); + spline1->SetLineColor(kMagenta+1); + spline1->SetLineWidth(1); + //spline1->Draw("SAME"); + + bin_max_hist_y1 = hEmpty2->GetMaximumBin(); + bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1); + max_spline1 = spline1->Eval(bins_center_x1); + + bin_min_hist_y1 = hEmpty2->GetMinimumBin(); + bins_center_xmin1 = hEmpty2->GetBinCenter(bin_min_hist_y1); + min_spline1 = spline1->Eval(bins_center_xmin1); + + for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){ + y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000)); + if (y_spline1 > (max_spline1+min_spline1)/2){ + x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000); + break; + } + } + for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){ + y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000)); + if (y_spline_r1 < (max_spline1+min_spline1)/2){ + x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000); + break; + } + } + + fwhm_spline1 = x_spline_r1 - x_spline_l1; + fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1); + } + fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1))); + } + } + return fwhm_spline; + return fwhmError; +} diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h index dcd7683f0d1..27f1448513e 100644 --- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h +++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h @@ -193,6 +193,11 @@ class AliBalancePsi : public TObject { Double_t &sigma, Double_t &sigmaError, Double_t &skewness, Double_t &skewnessError, Double_t &kurtosis, Double_t &kurtosisError); + + //++++++++++++++++++// + Double_t GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist, + Double_t &fwhm_spline, Double_t &fwhmError); + //++++++++++++++++++// TH2D *GetQAHistHBTbefore() {return fHistHBTbefore;} TH2D *GetQAHistHBTafter() {return fHistHBTafter;} diff --git a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C index d8a37294044..5613ff5449c 100644 --- a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C +++ b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C @@ -1,5 +1,5 @@ const Int_t numberOfCentralityBins = 12; -TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"}; +TString centralityArray[numberOfCentralityBins] = {"0-80","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"}; const Int_t gRebin = 1; @@ -933,7 +933,8 @@ void drawProjections(TH2D *gHistBalanceFunction2D = 0x0, TString eventClass = "Centrality", Bool_t bRootMoments = kFALSE, Bool_t kUseZYAM = kFALSE, - Bool_t bReduceRangeForMoments = kFALSE) { + Bool_t bReduceRangeForMoments = kFALSE, + Bool_t bFWHM = kFALSE) { //Macro that draws the balance functions PROJECTIONS //for each centrality bin for the different pT of trigger and //associated particles @@ -1292,12 +1293,37 @@ void drawProjections(TH2D *gHistBalanceFunction2D = 0x0, fileWeightedMean << " " << weightedMean << " " <SetFillColor(10); c2->SetHighLightColor(10); c2->SetLeftMargin(0.15); gHistBalanceFunctionSubtracted->DrawCopy("E"); + + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++// + Double_t fwhm_spline = 0.; + Double_t fwhmError = 0.; + if (bFWHM){ + AliBalancePsi *bHelper_1 = new AliBalancePsi; + bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError); + Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError); + + //store and in txt files FWHM + TString sigmaFileName = filename; + if(kProjectInEta) { + sigmaFileName = "deltaEtaProjection_FWHM.txt"; + } + else{ + sigmaFileName = "deltaPhiProjection_FWHM.txt"; + } + ofstream fileSigma(sigmaFileName.Data(),ios::app); + fileSigma << " " << fwhm_spline << " " <SetLineColor(2); + gHistBalanceFunctionSubtracted->SetMarkerColor(2); + gHistBalanceFunctionSubtracted->Draw("SAME"); + } + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++// TLatex *latex = new TLatex(); latex->SetNDC(); @@ -1340,8 +1366,8 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h", TString eventClass = "Multiplicity", Bool_t bRootMoments = kFALSE, Bool_t bFullPhiForEtaProjection = kFALSE, - Bool_t bReduceRangeForMoments = kFALSE -) { + Bool_t bReduceRangeForMoments = kFALSE, + Bool_t bFWHM = kFALSE) { //Macro that draws the BF distributions for each centrality bin //for reaction plane dependent analysis //Author: Panos.Christakoglou@nikhef.nl @@ -1479,10 +1505,11 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h", eventClass.Data(), bRootMoments, kFALSE, - bReduceRangeForMoments); + bReduceRangeForMoments, + bFWHM); } else{ - drawProjections(gHistBalanceFunction2D, + drawProjections(gHistBalanceFunction2D, kTRUE, 1,36, gCentrality, @@ -1491,9 +1518,10 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h", ptAssociatedMin,ptAssociatedMax, kTRUE, eventClass.Data(), - bRootMoments, - kFALSE, - bReduceRangeForMoments); + bRootMoments, + kFALSE, + bReduceRangeForMoments, + bFWHM); } drawProjections(gHistBalanceFunction2D, @@ -1507,7 +1535,8 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h", eventClass.Data(), bRootMoments, kFALSE, - bReduceRangeForMoments); + bReduceRangeForMoments, + bFWHM); TString outFileName = filename; outFileName.ReplaceAll("correlationFunction","balanceFunction2D"); @@ -1643,17 +1672,18 @@ sFileName[iDir] += momDirectory; hBFOut->Scale(1./entriesOut); drawProjections(hBFOut, - kTRUE, - 1,36, - gCentrality, - psiMin,psiMax, - ptTriggerMin,ptTriggerMax, - ptAssociatedMin,ptAssociatedMax, - kTRUE, - eventClass, - kTRUE, - kFALSE, - bReduceRangeForMoments); + kTRUE, + 1,36, + gCentrality, + psiMin,psiMax, + ptTriggerMin,ptTriggerMax, + ptAssociatedMin,ptAssociatedMax, + kTRUE, + eventClass, + kTRUE, + kFALSE, + bReduceRangeForMoments, + bFWHM); drawProjections(hBFOut, kFALSE, @@ -1666,7 +1696,8 @@ sFileName[iDir] += momDirectory; eventClass.Data(), kTRUE, kFALSE, - bReduceRangeForMoments); + bReduceRangeForMoments, + bFWHM); // output TString outFileName = "balanceFunction2D."; -- 2.39.3