#include "TGraphErrors.h"\r
#include "TH1F.h"\r
#include "TH2F.h"\r
+#include "TH3F.h" \r
#include "TH2D.h" \r
#include "TH3D.h"\r
#include "TArrayF.h"\r
fHistEta(0),\r
fHistRapidity(0),\r
fHistPhi(0),\r
+ fHistEtaPhiPos(0), \r
+ fHistEtaPhiNeg(0), \r
fHistPhiBefore(0),\r
fHistPhiAfter(0),\r
fHistPhiPos(0),\r
// delete fHistPt;\r
// delete fHistEta;\r
// delete fHistPhi;\r
+ // delete fHistEtaPhiPos; \r
+ // delete fHistEtaPhiNeg;\r
// delete fHistV0M;\r
}\r
\r
fList->Add(fHistRapidity);\r
fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhi);\r
+ fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
+ fList->Add(fHistEtaPhiPos); \r
+ fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
+ fList->Add(fHistEtaPhiNeg);\r
fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
fList->Add(fHistPhiBefore);\r
fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
\r
//=======================================correction\r
Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
fHistPt->Fill(vPt,gCentrality);\r
fHistEta->Fill(vEta,gCentrality);\r
fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
fHistRapidity->Fill(vY,gCentrality);\r
if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
fHistPt->Fill(vPt,gCentrality);\r
fHistEta->Fill(vEta,gCentrality);\r
fHistPhi->Fill(vPhi,gCentrality);\r
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
//fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
fHistRapidity->Fill(vY,gCentrality);\r
//if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
class TList;\r
class TH1F;\r
class TH2F;\r
+class TH3F; \r
class TF1;\r
class TH3D;\r
\r
TH2F *fHistEta;//pseudorapidity (QA histogram)\r
TH2F *fHistRapidity;//rapidity (QA histogram)\r
TH2F *fHistPhi;//phi (QA histogram)\r
+ TH3F *fHistEtaPhiPos;//eta-phi pos particles (QA histogram) \r
+ TH3F *fHistEtaPhiNeg;//eta-phi neg particles (QA histogram)\r
TH2F *fHistPhiBefore;//phi before v2 afterburner (QA histogram)\r
TH2F *fHistPhiAfter;//phi after v2 afterburner (QA histogram)\r
TH2F *fHistPhiPos;//phi for positive particles (QA histogram)\r
return gHistBalanceFunctionHistogram;
}
+//____________________________________________________________________//
+TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
+ Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax,
+ AliBalancePsi *bfMix) {
+ //Returns the BF histogram in Delta eta OR Delta phi,
+ //after dividing each correlation function by the Event Mixing one
+ // (But the division is done here in 2D, this was basically done to check the Integral)
+ //extracted from the 6 AliTHn objects
+ //(private members) of the AliBalancePsi class.
+ //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
+ //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
+ TString histName = "gHistBalanceFunctionHistogram1D";
+
+ if(!bfMix){
+ AliError("balance function object for event mixing not available");
+ return NULL;
+ }
+
+ // Psi_2
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+
+
+ // ============================================================================================
+ // the same for event mixing
+ AliTHn *fHistPMix = bfMix->GetHistNp();
+ AliTHn *fHistNMix = bfMix->GetHistNn();
+ AliTHn *fHistPNMix = bfMix->GetHistNpn();
+ AliTHn *fHistNPMix = bfMix->GetHistNnp();
+ AliTHn *fHistPPMix = bfMix->GetHistNpp();
+ AliTHn *fHistNNMix = bfMix->GetHistNnn();
+
+ // Psi_2
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+ // ============================================================================================
+
+
+ //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
+
+ // Project into the wanted space (1st: analysis step, 2nd: axis)
+ TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
+ TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
+ TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
+ TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
+ TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
+ TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
+
+ // ============================================================================================
+ // the same for event mixing
+ TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
+ TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
+ TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
+ TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
+ TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
+ TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
+ // ============================================================================================
+
+ TH1D *gHistBalanceFunctionHistogram = 0x0;
+ if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
+
+ hTemp1->Sumw2();
+ hTemp2->Sumw2();
+ hTemp3->Sumw2();
+ hTemp4->Sumw2();
+ hTemp1Mix->Sumw2();
+ hTemp2Mix->Sumw2();
+ hTemp3Mix->Sumw2();
+ hTemp4Mix->Sumw2();
+
+ hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp3->Scale(1./hTemp5->GetEntries());
+ hTemp2->Scale(1./hTemp6->GetEntries());
+ hTemp4->Scale(1./hTemp6->GetEntries());
+
+ hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
+ hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
+
+ hTemp1->Divide(hTemp1Mix);
+ hTemp2->Divide(hTemp2Mix);
+ hTemp3->Divide(hTemp3Mix);
+ hTemp4->Divide(hTemp4Mix);
+
+ // now only project on one axis
+ TH1D *h1DTemp1 = NULL;
+ TH1D *h1DTemp2 = NULL;
+ TH1D *h1DTemp3 = NULL;
+ TH1D *h1DTemp4 = NULL;
+ if(!bPhi){
+ h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
+ h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
+ h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
+ h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
+ }
+ else{
+ h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
+ h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
+ h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
+ h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
+ }
+
+ gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
+ gHistBalanceFunctionHistogram->Reset();
+ if(!bPhi){
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
+ }
+ else{
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
+ }
+
+ h1DTemp1->Add(h1DTemp3,-1.);
+ h1DTemp2->Add(h1DTemp4,-1.);
+
+ gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
+ gHistBalanceFunctionHistogram->Scale(0.5);
+ }
+
+ return gHistBalanceFunctionHistogram;
+}
+
//____________________________________________________________________//
TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
return gHistNN;
}
+//____________________________________________________________________//
+
+Bool_t AliBalancePsi::GetMomentsAnalytical(TH1D* gHist,
+ Double_t &mean, Double_t &meanError,
+ Double_t &sigma, Double_t &sigmaError,
+ Double_t &skewness, Double_t &skewnessError,
+ Double_t &kurtosis, Double_t &kurtosisError) {
+ //
+ // helper method to calculate the moments and errors of a TH1D anlytically
+ //
+
+ Bool_t success = kFALSE;
+ mean = -1.;
+ meanError = -1.;
+ sigma = -1.;
+ sigmaError = -1.;
+ skewness = -1.;
+ skewnessError = -1.;
+ kurtosis = -1.;
+ kurtosisError = -1.;
+
+ if(gHist){
+
+ // ----------------------------------------------------------------------
+ // basic parameters of histogram
+
+ Int_t fNumberOfBins = gHist->GetNbinsX();
+ // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
+
+
+ // ----------------------------------------------------------------------
+ // first calculate the mean
+
+ Double_t fWeightedAverage = 0.;
+ Double_t fNormalization = 0.;
+
+ for(Int_t i = 1; i <= fNumberOfBins; i++) {
+
+ fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
+ fNormalization += gHist->GetBinContent(i);
+
+ }
+
+ mean = fWeightedAverage / fNormalization;
+
+
+ // ----------------------------------------------------------------------
+ // then calculate the higher moments
+
+ Double_t fDelta = 0.;
+ Double_t fDelta2 = 0.;
+ Double_t fDelta3 = 0.;
+ Double_t fDelta4 = 0.;
+
+ for(Int_t i = 1; i <= fNumberOfBins; i++) {
+ fDelta += gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean;
+ fDelta2 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),2);
+ fDelta3 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),3);
+ fDelta4 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),4);
+ }
+
+ sigma = fDelta2 / fNormalization;
+ skewness = fDelta3 / fNormalization / TMath::Power(sigma,3);
+ kurtosis = fDelta4 / fNormalization / TMath::Power(sigma,4) - 3;
+
+ success = kTRUE;
+ }
+
+
+ return success;
+}
+
+
//____________________________________________________________________//
Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) {
//
Double_t ptAssociatedMin=-1.,
Double_t ptAssociatedMax=-1.,
AliBalancePsi *bfMix=NULL);
+
+ TH1D *GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
+ Double_t psiMin, Double_t psiMax,
+ Double_t ptTriggerMin=-1.,
+ Double_t ptTriggerMax=-1.,
+ Double_t ptAssociatedMin=-1.,
+ Double_t ptAssociatedMax=-1.,
+ AliBalancePsi *bfMix=NULL);
+
+ Bool_t GetMomentsAnalytical(TH1D* gHist,
+ Double_t &mean, Double_t &meanError,
+ Double_t &sigma, Double_t &sigmaError,
+ Double_t &skewness, Double_t &skewnessError,
+ Double_t &kurtosis, Double_t &kurtosisError);
TH2D *GetQAHistHBTbefore() {return fHistHBTbefore;}
TH2D *GetQAHistHBTafter() {return fHistHBTafter;}
Double_t ptAssociatedMax = -1.,
Bool_t k2pMethod = kFALSE,
Bool_t k2pMethod2D = kFALSE,
- TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
+ TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
+ Bool_t bRootMoments = kTRUE)
{
//Macro that draws the BF distributions for each centrality bin
//for reaction plane dependent analysis
gCentralityEstimator,
psiMin,psiMax,
ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
- k2pMethod,k2pMethod2D,eventClass);
+ k2pMethod,k2pMethod2D,eventClass,bRootMoments);
}
//______________________________________________________//
Double_t psiMin, Double_t psiMax,
Double_t ptTriggerMin, Double_t ptTriggerMax,
Double_t ptAssociatedMin, Double_t ptAssociatedMax,
- Bool_t k2pMethod = kFALSE,Bool_t k2pMethod2D = kFALSE, TString eventClass="EventPlane") {
+ Bool_t k2pMethod = kFALSE,Bool_t k2pMethod2D = kFALSE, TString eventClass="EventPlane",Bool_t bRootMoments=kTRUE) {
gROOT->LoadMacro("~/SetPlotStyle.C");
SetPlotStyle();
gStyle->SetPalette(1,0);
}
//Raw balance function
- if(k2pMethod){
+ if(k2pMethod && !k2pMethod2D){
if(bMixed){
gHistBalanceFunction = b->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
}
return;
}
}
- else if(k2pMethod2D){
+ else if(k2pMethod && k2pMethod2D){
if(bMixed){
if(gDeltaEtaDeltaPhi==1) //Delta eta
gHistBalanceFunction = b->GetBalanceFunction1DFrom2D2pMethod(0,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
//gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
//Mixed balance function
- if(k2pMethod){
+ if(k2pMethod && !k2pMethod2D){
if(bMixed)
gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
else{
return;
}
}
- else if(k2pMethod2D){
+ else if(k2pMethod && k2pMethod2D){
if(bMixed){
if(gDeltaEtaDeltaPhi==1) //Delta eta
gHistBalanceFunctionMixed = bMixed->GetBalanceFunction1DFrom2D2pMethod(0,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
//GetWeightedMean(gHistBalanceFunctionShuffled);
TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
- meanLatex = "#mu = ";
- meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
- meanLatex += " #pm ";
- meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
-
- rmsLatex = "#sigma = ";
- rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
- rmsLatex += " #pm ";
- rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
-
- skewnessLatex = "S = ";
- skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
- skewnessLatex += " #pm ";
- skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
-
- kurtosisLatex = "K = ";
- kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
- kurtosisLatex += " #pm ";
- kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
- Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
- Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
- Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
- Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
+
+ if(bRootMoments){
+ meanLatex = "#mu = ";
+ meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
+ meanLatex += " #pm ";
+ meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
+
+ rmsLatex = "#sigma = ";
+ rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
+ rmsLatex += " #pm ";
+ rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
+
+ skewnessLatex = "S = ";
+ skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
+ skewnessLatex += " #pm ";
+ skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
+
+ kurtosisLatex = "K = ";
+ kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
+ kurtosisLatex += " #pm ";
+ kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
+ Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
+ Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
+ Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
+ Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
+ }
+ // calculate the moments by hand
+ else{
+
+ Double_t meanAnalytical, meanAnalyticalError;
+ Double_t sigmaAnalytical, sigmaAnalyticalError;
+ Double_t skewnessAnalytical, skewnessAnalyticalError;
+ Double_t kurtosisAnalytical, kurtosisAnalyticalError;
+
+ b->GetMomentsAnalytical(gHistBalanceFunctionSubtracted,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
+
+ meanLatex = "#mu = ";
+ meanLatex += Form("%.3f",meanAnalytical);
+ meanLatex += " #pm ";
+ meanLatex += Form("%.3f",meanAnalyticalError);
+
+ rmsLatex = "#sigma = ";
+ rmsLatex += Form("%.3f",sigmaAnalytical);
+ rmsLatex += " #pm ";
+ rmsLatex += Form("%.3f",sigmaAnalyticalError);
+
+ skewnessLatex = "S = ";
+ skewnessLatex += Form("%.3f",skewnessAnalytical);
+ skewnessLatex += " #pm ";
+ skewnessLatex += Form("%.3f",skewnessAnalyticalError);
+
+ kurtosisLatex = "K = ";
+ kurtosisLatex += Form("%.3f",kurtosisAnalytical);
+ kurtosisLatex += " #pm ";
+ kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
+ Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
+ Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
+ Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
+ Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
+ }
TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
c2->SetFillColor(10);