]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
reverting changes that were made by accident in version 60358 AND possibilty to calcu...
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Mar 2013 18:17:22 +0000 (18:17 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Mar 2013 18:17:22 +0000 (18:17 +0000)
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
PWGCF/EBYE/macros/drawBalanceFunctionPsi.C

index 777ef08546f9fb946186699836a06b7aa6717238..25ac33adfd4e349fb59f69ea2d53e7cf10bb0cdc 100755 (executable)
@@ -5,6 +5,7 @@
 #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
@@ -79,6 +80,8 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   fHistEta(0),\r
   fHistRapidity(0),\r
   fHistPhi(0),\r
+  fHistEtaPhiPos(0),            \r
+  fHistEtaPhiNeg(0), \r
   fHistPhiBefore(0),\r
   fHistPhiAfter(0),\r
   fHistPhiPos(0),\r
@@ -192,6 +195,8 @@ AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
   // delete fHistPt;\r
   // delete fHistEta;\r
   // delete fHistPhi;\r
+  // delete fHistEtaPhiPos;                     \r
+  // delete fHistEtaPhiNeg;\r
   // delete fHistV0M;\r
 }\r
 \r
@@ -312,6 +317,10 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   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
@@ -1046,6 +1055,8 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gC
       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
@@ -1226,6 +1237,8 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gC
       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
@@ -1319,6 +1332,8 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gC
       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
index 24f72e99ba736ffa7b91d743f8ee28f4a074d59b..c0322f031adf44719e89a3ccb5cd4f5d2ab51937 100755 (executable)
@@ -7,6 +7,7 @@
 class TList;\r
 class TH1F;\r
 class TH2F;\r
+class TH3F; \r
 class TF1;\r
 class TH3D;\r
 \r
@@ -190,6 +191,8 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {
   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
index 49ab29f7d6b29715d1fd5c0e7f4d1412e7de175f..68f6d565f36760863ce6d38a28d97c0003930282 100644 (file)
@@ -986,6 +986,179 @@ TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
   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, 
@@ -1237,6 +1410,79 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(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) { 
   //
index 0891b66e2df82279b71ff79319d6cb2e1ed2e1d0..2b5ec67895c332454e43cfe29f365f76df19e3d4 100644 (file)
@@ -145,6 +145,20 @@ class AliBalancePsi : public TObject {
                                                   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;}
index fe55e1052ab9ee2ddc994dd188c4fbc45e7071c9..fb55c4b8d880a0707d01541b94eebcd909a22728 100644 (file)
@@ -13,7 +13,8 @@ void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root",
                            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
@@ -46,7 +47,7 @@ void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root",
         gCentralityEstimator,
         psiMin,psiMax,
         ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
-        k2pMethod,k2pMethod2D,eventClass);  
+        k2pMethod,k2pMethod2D,eventClass,bRootMoments);  
 }
 
 //______________________________________________________//
@@ -206,7 +207,7 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
          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);
@@ -384,7 +385,7 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
   }
 
   //Raw balance function
-  if(k2pMethod){ 
+  if(k2pMethod && !k2pMethod2D){ 
     if(bMixed){
       gHistBalanceFunction = b->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
     }
@@ -393,7 +394,7 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
       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);
@@ -439,7 +440,7 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
   //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{
@@ -447,7 +448,7 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
       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);
@@ -531,29 +532,66 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
   //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);