#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;
}
}
binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
}
- Double_t binPsiLowEdge = 0.;
- Double_t binPsiUpEdge = 0.;
- Double_t binVertexLowEdge = 0.;
- Double_t binVertexUpEdge = 0.;
TH2D* h1 = NULL;
TH2D* h2 = NULL;
cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
- // determine the bin edges for this bin
- binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
- binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
- if(fVertexBinning){
- binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
- binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
- }
-
// Psi_2
fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
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 ptAssociatedMin,
Double_t ptAssociatedMax,
AliBalancePsi *bMixed,
- Bool_t normToTrig) {
+ Bool_t normToTrig,
+ Double_t normalizationRangePhi) {
// Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
// does the division by event mixing inside,
TH2D *fMixed = NULL;
// ranges (in bins) for vertexZ and centrality (psi)
- Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
+ 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);
+ 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 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++){
- cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
+ //cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
// determine the bin edges for this bin
- binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
- binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
+ binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
+ binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
if(fVertexBinning){
- binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
- binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
+ binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
+ binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
}
// get the 2D histograms for the correct type
fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
}
-
+
if(fMixed && normToTrig && fMixed->Integral()>0){
// normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
- // do it only on away-side (due to two-track cuts)
- Double_t mixedNorm = fMixed->Integral(fMixed->GetXaxis()->FindBin(0-10e-5),fMixed->GetXaxis()->FindBin(0+10e-5),fMixed->GetNbinsY()/2+1,fMixed->GetNbinsY());
- mixedNorm /= 0.5 * fMixed->GetNbinsY() *(fMixed->GetXaxis()->FindBin(0.01) - fMixed->GetXaxis()->FindBin(-0.01) + 1);
+ // do it only on away-side (due to two-track cuts): pi +- pi/6.
+ Int_t binXmin = fMixed->GetXaxis()->FindBin(0-10e-5);
+ Int_t binXmax = fMixed->GetXaxis()->FindBin(0+10e-5);
+ Double_t binsX = (Double_t)(binXmax - binXmin + 1);
+ Int_t binYmin = fMixed->GetYaxis()->FindBin(TMath::Pi() - normalizationRangePhi);
+ Int_t binYmax = fMixed->GetYaxis()->FindBin(TMath::Pi() + normalizationRangePhi - 0.00001);
+ Double_t binsY = (Double_t)(binYmax - binYmin + 1);
+
+ Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
+ mixedNorm /= binsX * binsY;
// finite bin correction
Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
if(fSame && fMixed){
// then get the correlation function (divide fSame/fmixed)
- fSame->Divide(fMixed);
-
- // NEW averaging:
+
+ 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;
if(type=="PN" || type=="PP")
NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
else if(type=="NP" || type=="NN")
NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
+ else if(type=="ALL")
+ 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;
}
+
}
}
}
if(gHist){
- // OLD averaging:
- // average over number of bins nbinsVertex * nbinsPsi
- // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
-
- // NEW averaging:
- // average over number of triggers in each sub-bin
+ // averaging with number of triggers:
// first set to full range and then obtain number of all triggers
Double_t NTrigAll = 0;
if(type=="PN" || type=="PP"){
fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
}
+ 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);
+ 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;
}
// pt associated
- if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)){
fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
-
+ }
+
//0:step, 1: Delta eta, 2: Delta phi
TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
if(!gHistNN){
// divide by sum of + and - triggers
if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
- gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistN->Project(0,1)->Integral()));
+ gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral()));
//normalize to bin width
gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
}
//____________________________________________________________________//
-
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