Modifications for proper scaling to bin width, Event Mixing bins for centrality and...

diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
index 89e5f2b2492627607b21d9037f20385a052e3627..1a0aa7ef58f09b51fba9b2d47445af6933d3d309 100755 (executable)
--- a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
+++ b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
@@ -400,9 +400,14 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
     \r
     // centrality bins\r
-    Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
+    Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
     Double_t* centbins        = centralityBins;\r
     Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
+\r
+    // multiplicity bins\r
+    Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
+    Double_t* multbins        = multiplicityBins;\r
+    Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
     \r
     // Zvtx bins\r
     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
@@ -416,10 +421,20 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     \r
     // run the event mixing also in bins of event plane (statistics!)\r
     if(fRunMixingEventPlane){\r
-      fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
+      if(fEventClass=="Multiplicity"){\r
+       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
+      }\r
+      else{\r
+       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
+      }\r
     }\r
     else{\r
-      fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
+      if(fEventClass=="Multiplicity"){\r
+       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
+      }\r
+      else{\r
+       fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
+      }\r
     }\r
   }\r
 \r
@@ -590,17 +605,17 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
 \r
   // get the reaction plane\r
   gReactionPlane = GetEventPlane(eventMain);\r
-  fHistEventPlane->Fill(gReactionPlane,gCentrality);\r
+  fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
   if(gReactionPlane < 0){\r
     return;\r
   }\r
   \r
   // get the accepted tracks in main event\r
-  TObjArray *tracksMain = GetAcceptedTracks(eventMain,gCentrality,gReactionPlane);\r
+  TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 \r
   //multiplicity cut (used in pp)\r
-  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,gCentrality);\r
+  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
   if(fUseMultiplicity) {\r
     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
       return;\r
@@ -609,7 +624,7 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
   // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
   TObjArray* tracksShuffled = NULL;\r
   if(fRunShuffling){\r
-    tracksShuffled = GetShuffledTracks(tracksMain,gCentrality);\r
+    tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
   }\r
   \r
   // Event mixing \r
@@ -631,10 +646,10 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
       //    of 1./nMix can be applied.\r
       \r
-      AliEventPool* pool = fPoolMgr->GetEventPool(gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
+      AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
       \r
       if (!pool){\r
-       AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", gCentrality, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
+       AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
       }\r
       else{\r
        \r

diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
index 23fdaab60f3f9cc1d37ec256cd8a954377ba449b..b4ee4424b051159a99c1e46a531c2808422772ea 100644 (file)
--- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
+++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
@@ -177,16 +177,16 @@ void AliBalancePsi::InitHistograms() {
   ***********************************************************/
    
   //--- Multiplicity Bins ------------------------------------
-    const Int_t kMultBins = 8;
+    const Int_t kMultBins = 10;
     //A first rough attempt at four bins
-    Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80};
+    Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80,100,100000};
   //----------------------------------------------------------
     
   //--- Centrality Bins --------------------------------------
     const Int_t kNCentralityBins       = 9;
-    const Int_t kNCentralityBinsVertex = 26;
+    const Int_t kNCentralityBinsVertex = 15;
     Double_t centralityBins[kNCentralityBins+1]             = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
-    Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.};
+    Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.};
   //----------------------------------------------------------
     
   //--- Event Plane Bins -------------------------------------
@@ -266,9 +266,9 @@ void AliBalancePsi::InitHistograms() {
 
   // pt(trigger-associated)
   const Int_t kNPtBins       = 16;
-  const Int_t kNPtBinsVertex = 5;
+  const Int_t kNPtBinsVertex = 6;
   Double_t ptBins[kNPtBins+1]             = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
-  Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0};
+  Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0,15.0};
 
   // coarse binning in case of vertex Z binning
   if(fVertexBinning){
@@ -796,7 +796,7 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
     gHistBalanceFunctionHistogram->Scale(0.5);
 
     //normalize to bin width
-    gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
+    gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
   }
 
   return gHistBalanceFunctionHistogram;
@@ -1044,7 +1044,7 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
     gHistBalanceFunctionHistogram->Scale(0.5);
 
     //normalize to bin width
-    gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
+    gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
   }
 
   return gHistBalanceFunctionHistogram;
@@ -1661,6 +1661,9 @@ TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
 
     gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
     gHistBalanceFunctionHistogram->Scale(0.5);
+
+    //normalize to bin width
+    gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
   }
 
   return gHistBalanceFunctionHistogram;
@@ -2258,9 +2261,9 @@ Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist,
     // kurtosisError = TMath::Sqrt(24./(normError));
 
     // use delta theorem paper (Luo - arXiv:1109.0593v1)
-    Double_t Lambda11 = (fMu4-1)*sigma*sigma/(4*normError);
-    Double_t Lambda22 = (9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError;
-    Double_t Lambda33 = (-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError;
+    Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
+    Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
+    Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
     //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
     //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
     //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
@@ -2268,11 +2271,14 @@ Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist,
     // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
     // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
 
-    meanError        = sigma / TMath::Sqrt(normError); 
+    if (TMath::Sqrt(normError) != 0){
+      meanError        = sigma / TMath::Sqrt(normError); 
+    }
+    else return -999;
     sigmaError       = TMath::Sqrt(Lambda11);
     skewnessError    = TMath::Sqrt(Lambda22);
     kurtosisError    = TMath::Sqrt(Lambda33);
-
+    
     
     success = kTRUE;    
   }

diff --git a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
index deac0be3b7c9e16b84e9a24ef3f4d14086a651d2..a7aa0b9ac4561b6c7100106a18fa4445f9924b39 100644 (file)
--- 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-10","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-100","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;
@@ -530,7 +530,12 @@ void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
   newFileName += Form("%.1f",ptAssociatedMax); 
   if(k2pMethod) newFileName += "_2pMethod";
-  newFileName += ".root";
+
+  // newFileName += "_";
+  // newFileName += Form("%.1f",psiMin); 
+  // newFileName += "-"; 
+  // newFileName += Form("%.1f",psiMax); 
+  // newFileName += ".root";
 
   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
   fOutput->cd();
@@ -650,7 +655,13 @@ void fitbalanceFunction(Int_t gCentrality = 1,
   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
   newFileName += Form("%.1f",ptAssociatedMax); 
   if(k2pMethod) newFileName += "_2pMethod";
+
+  newFileName += "_";
+  newFileName += Form("%.1f",psiMin); 
+  newFileName += "-"; 
+  newFileName += Form("%.1f",psiMax); 
   newFileName += ".root";
+
   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
   gHist->Write();
   gHistResidual->Write();
@@ -699,6 +710,11 @@ void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
   filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
   filename += Form("%.1f",ptAssociatedMax);  
   if(k2pMethod) filename += "_2pMethod";
+
+  filename += "_";
+  filename += Form("%.1f",psiMin); 
+  filename += "-"; 
+  filename += Form("%.1f",psiMax); 
   filename += ".root";  
 
   //Open the file
@@ -918,6 +934,16 @@ void drawProjections(const char* lhcPeriod = "LHC10h",
   TGaxis::SetMaxDigits(3);
 
   //first we need some libraries
+  gSystem->Load("libTree");
+  gSystem->Load("libGeom");
+  gSystem->Load("libVMC");
+  gSystem->Load("libXMLIO");
+  gSystem->Load("libPhysics");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libESD");
+  gSystem->Load("libAOD");
+
   gSystem->Load("libANALYSIS.so");
   gSystem->Load("libANALYSISalice.so");
   gSystem->Load("libEventMixing.so");
@@ -949,6 +975,11 @@ void drawProjections(const char* lhcPeriod = "LHC10h",
   filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
   filename += Form("%.1f",ptAssociatedMax); 
   if(k2pMethod) filename += "_2pMethod";
+
+  // filename += "_";
+  // filename += Form("%.1f",psiMin); 
+  // filename += "-"; 
+  // filename += Form("%.1f",psiMax);
   filename += ".root";
 
   //Open the file
@@ -998,15 +1029,22 @@ void drawProjections(const char* lhcPeriod = "LHC10h",
   TH1D *gHistBalanceFunctionSubtracted = NULL;
   TH1D *gHistBalanceFunctionMixed      = NULL;
 
+  TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
+  TH1D *gHistBalanceFunctionMixed_scale      = NULL;
+
   if(kProjectInEta){
     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX());
+    gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX());
+    gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");  
   }
   else{
     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY());
+    gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY());
+    gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");  
   }
@@ -1023,7 +1061,7 @@ void drawProjections(const char* lhcPeriod = "LHC10h",
   c1->SetHighLightColor(10);
   c1->SetLeftMargin(0.15);
   gHistBalanceFunctionSubtracted->DrawCopy("E");
-  gHistBalanceFunctionMixed->DrawCopy("E, SAME");
+  gHistBalanceFunctionMixed->DrawCopy("ESAME");
   
   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
   legend->SetTextSize(0.045); 
@@ -1064,6 +1102,10 @@ void drawProjections(const char* lhcPeriod = "LHC10h",
   pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
   pngName += Form("%.1f",ptAssociatedMax); 
   if(k2pMethod) pngName += "_2pMethod";
+  
+  pngName += "_"; 
+  pngName += Form("%.1f",psiMin); pngName += "-"; 
+  pngName += Form("%.1f",psiMax);
   pngName += ".png";
 
   c1->SaveAs(pngName.Data());