#include <TObjArray.h>
#include <TGraphErrors.h>
#include <TString.h>
+#include <TSpline.h>
+#include <TRandom3.h>
#include "AliVParticle.h"
#include "AliMCParticle.h"
}//resonance cut
// HBT like cut
- if(fHBTCut){ // VERSION 3 (all pairs)
- //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
+ //if(fHBTCut){ // VERSION 3 (all pairs)
+ if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
//if( dphi < 3 || deta < 0.01 ){ // VERSION 1
// continue;
const Float_t kLimit = fHBTCutValue * 3;
Float_t dphistarminabs = 1e5;
- Float_t dphistarmin = 1e5;
+ //Float_t dphistarmin = 1e5;
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 dphistarabs = TMath::Abs(dphistar);
if (dphistarabs < dphistarminabs) {
- dphistarmin = dphistar;
+ //dphistarmin = dphistar;
dphistarminabs = dphistarabs;
}
}
return gHistBalanceFunctionHistogram;
}
+//____________________________________________________________________//
+TH1D *AliBalancePsi::GetTriggers(TString type,
+ Double_t psiMin,
+ Double_t psiMax,
+ Double_t vertexZMin,
+ Double_t vertexZMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax
+ ) {
+
+ // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
+ // does the division by event mixing inside,
+ // and averaging over several vertexZ and centrality bins
+
+ // security checks
+ if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
+ AliError("Only types allowed: PN,NP,PP,NN,ALL");
+ return NULL;
+ }
+
+ TH1D *gHist = NULL;
+
+ // ranges (in bins) for vertexZ and centrality (psi)
+ Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
+ Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
+ Int_t binVertexMin = 0;
+ Int_t binVertexMax = 0;
+ if(fVertexBinning){
+ binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
+ binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
+ }
+ Double_t binPsiLowEdge = 0.;
+ Double_t binPsiUpEdge = 1000.;
+ Double_t binVertexLowEdge = 0.;
+ Double_t binVertexUpEdge = 1000.;
+
+
+ // first set to range and then obtain histogram
+ if(type=="PN" || type=="PP"){
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistP->Project(0,1);
+ }
+ else if(type=="NP" || type=="NN"){
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistN->Project(0,1);
+ }
+ else if(type=="ALL"){
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistN->Project(0,1);
+ gHist->Add((TH1D*)fHistP->Project(0,1));
+ }
+
+ return gHist;
+}
//____________________________________________________________________//
TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
Double_t binVertexLowEdge = 0.;
Double_t binVertexUpEdge = 1000.;
+ // if event mixing is empty, we can not add that sub bin
+ // need to record the number of triggers for correct normalization
+ Double_t nTrigSubBinEmpty = 0.;
+
// loop over all bins
for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
if(fSame && fMixed){
// then get the correlation function (divide fSame/fmixed)
- fSame->Divide(fMixed);
-
+
+ if(fMixed->Integral()>0)
+ fSame->Divide(fMixed);
+
// averaging with number of triggers:
// average over number of triggers in each sub-bin
Double_t NTrigSubBin = 0;
NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
fSame->Scale(NTrigSubBin);
- // for the first: clone
- if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
- gHist = (TH2D*)fSame->Clone();
+ // only if event mixing has enough statistics
+ if(fMixed->Integral()>0){
+
+ // for the first: clone
+ if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
+ gHist = (TH2D*)fSame->Clone();
+ }
+ else{ // otherwise: add for averaging
+ gHist->Add(fSame);
+ }
}
- else{ // otherwise: add for averaging
- gHist->Add(fSame);
+
+ // otherwise record the number of triggers for correct normalization
+ else{
+ nTrigSubBinEmpty += NTrigSubBin;
}
+
}
}
}
fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
}
+
+ // subtract number of triggers with empty sub bins for correct normalization
+ NTrigAll -= nTrigSubBinEmpty;
+
gHist->Scale(1./NTrigAll);
-
}
return gHist;
}
//____________________________________________________________________//
-
Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
Double_t &mean, Double_t &meanError,
Double_t &sigma, Double_t &sigmaError,
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; f<repeat; f++){ //100
+ for (Int_t h=0; h<Nbin+1; h++){
+ dmeans = gHist->GetBinContent(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; f<repeat; f++){ //100
+ for (Int_t h=0; h<Nbin+1; h++){
+ dmeans = gHist->GetBinContent(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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff