From a38fd7f358bfbf7bc466479c475211b429876115 Mon Sep 17 00:00:00 2001 From: pchrist Date: Sun, 20 May 2012 11:30:52 +0000 Subject: [PATCH] Adding the pt trigger and associated bins in the Psi dependent bf analysis --- PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx | 160 +++++++-- PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h | 35 +- .../macros/AddTaskBalancePsiCentralityTrain.C | 42 +-- PWGCF/EBYE/macros/drawBalanceFunctionPsi.C | 323 ++++++++++++++++++ .../EBYE/macros/drawCorrelationFunctionPsi.C | 224 ++++++++++++ 5 files changed, 704 insertions(+), 80 deletions(-) create mode 100644 PWGCF/EBYE/macros/drawBalanceFunctionPsi.C create mode 100644 PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx index 9f744e112ce..1e6a3c8c02c 100644 --- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx +++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx @@ -129,10 +129,10 @@ void AliBalancePsi::InitHistograms() { axisTitlePair[1] = "#phi - #Psi_{2} (#circ)"; // eta - const Int_t kNEtaBins = 20; + const Int_t kNEtaBins = 10; Double_t etaBins[kNEtaBins+1]; for(Int_t i = 0; i < kNEtaBins+1; i++) - etaBins[i] = -1.0 + i * 0.1; + etaBins[i] = -1.0 + i * 0.2; iBinSingle[2] = kNEtaBins; dBinsSingle[2] = etaBins; axisTitleSingle[2] = "#eta"; @@ -143,10 +143,10 @@ void AliBalancePsi::InitHistograms() { axisTitlePair[2] = "#eta"; */ // delta eta - const Int_t kNDeltaEtaBins = 20; + const Int_t kNDeltaEtaBins = 80; Double_t deltaEtaBins[kNDeltaEtaBins+1]; for(Int_t i = 0; i < kNDeltaEtaBins+1; i++) - deltaEtaBins[i] = -2.0 + i * 0.2; + deltaEtaBins[i] = -2.0 + i * 0.05; iBinPair[2] = kNDeltaEtaBins; dBinsPair[2] = deltaEtaBins; axisTitlePair[2] = "#Delta #eta"; @@ -175,37 +175,15 @@ void AliBalancePsi::InitHistograms() { axisTitlePair[8] = "#Delta y"; */ // phi - const Int_t kNPhiBins = 36; + const Int_t kNPhiBins = 18; Double_t phiBins[kNPhiBins+1]; for(Int_t i = 0; i < kNPhiBins+1; i++){ - phiBins[i] = 0.0 + i * 10.; + phiBins[i] = 0.0 + i * 20.; } iBinSingle[3] = kNPhiBins; dBinsSingle[3] = phiBins; axisTitleSingle[3] = "#phi (#circ)"; - - /*iBinPair[4] = kNPhiBins; - dBinsPair[4] = phiBins; - axisTitlePair[4] = "#phi (#circ)"; */ - // pt(trigger) - /*const Int_t kNPtBins = 100; - Double_t ptBins[kNPtBins+1]; - for(Int_t i = 0; i < kNPtBins+1; i++){ - ptBins[i] = 0.0 + i * 0.2; - } - iBinSingle[4] = kNPtBins; - dBinsSingle[4] = ptBins; - axisTitleSingle[4] = "p_{t}^{trig.} (GeV/c)"; */ - - /*iBinPair[5] = kNPtBins; - dBinsPair[5] = ptBins; - axisTitlePair[5] = "p_{t}^{trig.} (GeV/c)"; - - iBinPair[6] = kNPtBins; - dBinsPair[6] = ptBins; - axisTitlePair[6] = "p_{t}^{assoc.} (GeV/c)"; */ - // delta phi const Int_t kNDeltaPhiBins = 72; Double_t deltaPhiBins[kNDeltaPhiBins+1]; @@ -216,6 +194,24 @@ void AliBalancePsi::InitHistograms() { dBinsPair[3] = deltaPhiBins; axisTitlePair[3] = "#Delta #phi (#circ)"; + // pt(trigger-associated) + const Int_t kNPtBins = 50; + Double_t ptBins[kNPtBins+1]; + for(Int_t i = 0; i < kNPtBins+1; i++){ + ptBins[i] = 0.0 + i * 0.5; + } + iBinSingle[4] = kNPtBins; + dBinsSingle[4] = ptBins; + axisTitleSingle[4] = "p_{t}^{trig.} (GeV/c)"; + + iBinPair[4] = kNPtBins; + dBinsPair[4] = ptBins; + axisTitlePair[4] = "p_{t}^{trig.} (GeV/c)"; + + iBinPair[5] = kNPtBins; + dBinsPair[5] = ptBins; + axisTitlePair[5] = "p_{t}^{assoc.} (GeV/c)"; + // qside /*const Int_t kNQSideBins = 200; Double_t qSideBins[kNQSideBins+1]; @@ -373,7 +369,7 @@ void AliBalancePsi::CalculateBalance(Float_t fCentrality, trackVarsSingle[2] = eta1; //trackVarsSingle[3] = rap1; trackVarsSingle[3] = phi1; - //trackVarsSingle[5] = pt1; + trackVarsSingle[4] = pt1; //Printf("Track(a) %d - phi-Psi: %lf",i+1,trackVarsSingle[1]); //fill single particle histograms @@ -443,11 +439,11 @@ void AliBalancePsi::CalculateBalance(Float_t fCentrality, //trackVarsPair[2] = eta1; //trackVarsPair[3] = rap1; //trackVarsPair[4] = phi1; - //trackVarsPair[5] = pt1; - //trackVarsPair[6] = pt2; trackVarsPair[2] = deta; //trackVarsPair[8] = dy; trackVarsPair[3] = dphi; + trackVarsPair[4] = pt1; + trackVarsPair[5] = pt2; //trackVarsPair[10] = qSide; //trackVarsPair[11] = qOut; //trackVarsPair[12] = qLong; @@ -494,6 +490,8 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle, fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + //Printf("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) TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); @@ -544,9 +542,9 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle, hTemp2->Sumw2(); hTemp3->Sumw2(); hTemp4->Sumw2(); - hTemp1->Add(hTemp3,-2.); + hTemp1->Add(hTemp3,-1.); hTemp1->Scale(1./hTemp5->GetEntries()); - hTemp2->Add(hTemp4,-2.); + hTemp2->Add(hTemp4,-1.); hTemp2->Scale(1./hTemp6->GetEntries()); gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.); gHistBalanceFunctionHistogram->Scale(0.5); @@ -554,3 +552,99 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle, return gHistBalanceFunctionHistogram; } + +//____________________________________________________________________// +TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t centrMin, + Double_t centrMax, + Double_t psiMin, + Double_t psiMax) { + //Returns the 2D correlation function for (+-) pairs + // centrality: axis 0 + fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + + // Psi_2: axis 1 + fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + fHistPN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + + //0:step, 2: Delta eta, 3: Delta phi + TH2D *gHist = dynamic_cast(fHistPN->Project(0,2,3)); + //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries())); + //Printf("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,2,3)->GetEntries())); + if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0) + gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries())); + + return gHist; +} + +//____________________________________________________________________// +TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t centrMin, + Double_t centrMax, + Double_t psiMin, + Double_t psiMax) { + //Returns the 2D correlation function for (+-) pairs + // centrality: axis 0 + fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + + // Psi_2: axis 1 + fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + fHistNP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + + //0:step, 2: Delta eta, 3: Delta phi + TH2D *gHist = dynamic_cast(fHistNP->Project(0,2,3)); + //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries())); + //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries())); + if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0) + gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries())); + + return gHist; +} + +//____________________________________________________________________// +TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t centrMin, + Double_t centrMax, + Double_t psiMin, + Double_t psiMax) { + //Returns the 2D correlation function for (+-) pairs + // centrality: axis 0 + fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + + // Psi_2: axis 1 + fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + fHistPP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + + //0:step, 2: Delta eta, 3: Delta phi + TH2D *gHist = dynamic_cast(fHistPP->Project(0,2,3)); + //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries())); + //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries())); + if((Double_t)(fHistP->Project(0,2)->GetEntries())!=0) + gHist->Scale(1./(Double_t)(fHistP->Project(0,2)->GetEntries())); + + return gHist; +} + +//____________________________________________________________________// +TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t centrMin, + Double_t centrMax, + Double_t psiMin, + Double_t psiMax) { + //Returns the 2D correlation function for (+-) pairs + // centrality: axis 0 + fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(centrMin,centrMax); + + // Psi_2: axis 1 + fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + fHistNN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(psiMin,psiMax); + + //0:step, 2: Delta eta, 3: Delta phi + TH2D *gHist = dynamic_cast(fHistNN->Project(0,2,3)); + //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries())); + //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries())); + if((Double_t)(fHistN->Project(0,2)->GetEntries())!=0) + gHist->Scale(1./(Double_t)(fHistN->Project(0,2)->GetEntries())); + + return gHist; +} diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h index 10580c224d7..4a09c3b0a28 100644 --- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h +++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h @@ -15,18 +15,18 @@ #include #include "TString.h" +#include "AliTHn.h" + #define ANALYSIS_TYPES 7 #define MAXIMUM_NUMBER_OF_STEPS 1024 #define MAXIMUM_STEPS_IN_PSI 360 -class TGraphErrors; class TH1D; class TH2D; class TH3D; -class AliTHn; -const Int_t nTrackVariablesSingle = 4; // track variables in histogram (eta, phi, pTtrig, centrality) -const Int_t nTrackVariablesPair = 4; // track variables in histogram (dEta, dPhi, pT, pTtrig, centrality) +const Int_t nTrackVariablesSingle = 5; // track variables in histogram (centrality, phi-Psi2, eta, phi, pTtrig) +const Int_t nTrackVariablesPair = 6; // track variables in histogram (centrality, phi-Psi2, dEta, dPhi, pTtrig, ptAssociated) const TString gBFPsiAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"}; class AliBalancePsi : public TObject { @@ -60,6 +60,15 @@ class AliBalancePsi : public TObject { void CalculateBalance(Float_t fCentrality, Double_t gReactionPlane, vector **chargeVector); + TH2D *GetCorrelationFunctionPN(Double_t centrMin, Double_t centrMax, + Double_t psiMin, Double_t psiMax); + TH2D *GetCorrelationFunctionNP(Double_t centrMin, Double_t centrMax, + Double_t psiMin, Double_t psiMax); + TH2D *GetCorrelationFunctionPP(Double_t centrMin, Double_t centrMax, + Double_t psiMin, Double_t psiMax); + TH2D *GetCorrelationFunctionNN(Double_t centrMin, Double_t centrMax, + Double_t psiMin, Double_t psiMax); + AliTHn *GetHistNp() {return fHistP;} AliTHn *GetHistNn() {return fHistN;} AliTHn *GetHistNpn() {return fHistPN;} @@ -67,12 +76,18 @@ class AliBalancePsi : public TObject { AliTHn *GetHistNpp() {return fHistPP;} AliTHn *GetHistNnn() {return fHistNN;} - void SetHistNp(AliTHn *gHist) {fHistP = gHist;} - void SetHistNn(AliTHn *gHist) {fHistN = gHist;} - void SetHistNpn(AliTHn *gHist) {fHistPN = gHist;} - void SetHistNnp(AliTHn *gHist) {fHistNP = gHist;} - void SetHistNpp(AliTHn *gHist) {fHistPP = gHist;} - void SetHistNnn(AliTHn *gHist) {fHistNN = gHist;} + void SetHistNp(AliTHn *gHist) { + fHistP = gHist; }//fHistP->FillParent(); fHistP->DeleteContainers();} + void SetHistNn(AliTHn *gHist) { + fHistN = gHist; }//fHistN->FillParent(); fHistN->DeleteContainers();} + void SetHistNpn(AliTHn *gHist) { + fHistPN = gHist; }//fHistPN->FillParent(); fHistPN->DeleteContainers();} + void SetHistNnp(AliTHn *gHist) { + fHistNP = gHist; }//fHistNP->FillParent(); fHistNP->DeleteContainers();} + void SetHistNpp(AliTHn *gHist) { + fHistPP = gHist; }//fHistPP->FillParent(); fHistPP->DeleteContainers();} + void SetHistNnn(AliTHn *gHist) { + fHistNN = gHist; }//fHistNN->FillParent(); fHistNN->DeleteContainers();} TH1D *GetBalanceFunctionHistogram(Int_t iVariableSingle, Int_t iVariablePair, diff --git a/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C b/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C index 03b5530bff8..2b5ce5be92f 100644 --- a/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C +++ b/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C @@ -34,38 +34,6 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0., TString fileNameBase="AnalysisResults") { // Creates a balance function analysis task and adds it to the analysis manager. // Get the pointer to the existing analysis manager via the static access method. - TString centralityName(""); - centralityName+=Form("%.0f",centrMin); - centralityName+="-"; - centralityName+=Form("%.0f",centrMax); - centralityName+="_"; - centralityName+=Form("%s",centralityEstimator.Data()); - centralityName+="_"; - centralityName+=Form("vZ%.1f",vertexZ); - centralityName+="_"; - centralityName+=Form("DCAxy%.1f",DCAxy); - centralityName+="_"; - centralityName+=Form("DCAz%.1f",DCAz); - centralityName+="_Pt"; - centralityName+=Form("%.1f",ptMin); - centralityName+="-"; - centralityName+=Form("%.1f",ptMax); - centralityName+="_Eta"; - centralityName+=Form("%.1f",etaMin); - centralityName+="-"; - centralityName+=Form("%.1f",etaMax); - centralityName+="_Chi"; - centralityName+=Form("%.1f",maxTPCchi2); - centralityName+="_nClus"; - centralityName+=Form("%d",minNClustersTPC); - centralityName+="_Bit"; - centralityName+=Form("%d",AODfilterBit); - if(bCentralTrigger) centralityName+="_withCentralTrigger"; - - - - - TString outputFileName(fileNameBase); outputFileName.Append(".root"); @@ -165,11 +133,11 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0., // Get and connect other common input/output containers via the manager as below //============================================================================== TString outputFileName = AliAnalysisManager::GetCommonFileName(); - outputFileName += ":PWGCFEbyE.outputBalanceFunctionAnalysis"; - AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQA_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); - AliAnalysisDataContainer *coutBF = mgr->CreateContainer(Form("listBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); - if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer(Form("listBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); - if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer(Form("listQAPID_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); + outputFileName += ":PWGCFEbyE.outputBalanceFunctionPsiAnalysis"; + AliAnalysisDataContainer *coutQA = mgr->CreateContainer("listQAPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); + AliAnalysisDataContainer *coutBF = mgr->CreateContainer("listBFPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); + if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer("listBFPsiShuffled", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); + if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer("listQAPIDPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data()); mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer()); mgr->ConnectOutput(taskBF, 1, coutQA); diff --git a/PWGCF/EBYE/macros/drawBalanceFunctionPsi.C b/PWGCF/EBYE/macros/drawBalanceFunctionPsi.C new file mode 100644 index 00000000000..00c0bc4e041 --- /dev/null +++ b/PWGCF/EBYE/macros/drawBalanceFunctionPsi.C @@ -0,0 +1,323 @@ +const Int_t numberOfCentralityBins = 9; +TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"}; +Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.}; +Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.}; +TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"}; + +const Int_t gRebin = 1; +void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root", + Double_t psiMin = 0., Double_t psiMax = 7.5) { + //Macro that draws the BF distributions for each centrality bin + //for reaction plane dependent analysis + //Author: Panos.Christakoglou@nikhef.nl + //Load the PWG2ebye library + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + gSystem->Load("libEventMixing.so"); + gSystem->Load("libCORRFW.so"); + gSystem->Load("libPWGTools.so"); + gSystem->Load("libPWGCFebye.so"); + + //Prepare the objects and return them + TList *listBF = GetListOfObjects(filename); + TList *listBFShuffled = GetListOfObjects(filename,kTRUE); + if(!listBF) { + Printf("The TList object was not created"); + return; + } + else + draw(listBF,listBFShuffled,psiMin,psiMax); +} + +//______________________________________________________// +TList *GetListOfObjects(const char* filename, Bool_t kShuffling = kFALSE) { + //Get the TList objects (QA, bf, bf shuffled) + TList *listQA = 0x0; + TList *listBF = 0x0; + TList *listBFShuffling = 0x0; + + //Open the file + TFile *f = TFile::Open(filename); + if((!f)||(!f->IsOpen())) { + Printf("The file %s is not found. Aborting...",filename); + return listBF; + } + //f->ls(); + + TDirectoryFile *dir = dynamic_cast(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis")); + if(!dir) { + Printf("The TDirectoryFile is not found. Aborting...",filename); + return listBF; + } + //dir->ls(); + + TString listBFName; + if(!kShuffling) + listBFName = "listBF_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger"; + else if(kShuffling) + listBFName = "listBFShuffled_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger"; + listBF = dynamic_cast(dir->Get(listBFName.Data())); + //listBF->ls(); + + //Get the histograms + TString histoName; + if(!kShuffling) + histoName = "fHistPV0M"; + else if(kShuffling) + histoName = "fHistP_shuffleV0M"; + AliTHn *fHistP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistP) { + Printf("fHistP %s not found!!!",histoName.Data()); + break; + } + fHistP->FillParent(); fHistP->DeleteContainers(); + + if(!kShuffling) + histoName = "fHistNV0M"; + else if(kShuffling) + histoName = "fHistN_shuffleV0M"; + AliTHn *fHistN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistN) { + Printf("fHistN %s not found!!!",histoName.Data()); + break; + } + fHistN->FillParent(); fHistN->DeleteContainers(); + + if(!kShuffling) + histoName = "fHistPNV0M"; + else if(kShuffling) + histoName = "fHistPN_shuffleV0M"; + AliTHn *fHistPN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistPN) { + Printf("fHistPN %s not found!!!",histoName.Data()); + break; + } + fHistPN->FillParent(); fHistPN->DeleteContainers(); + + if(!kShuffling) + histoName = "fHistNPV0M"; + else if(kShuffling) + histoName = "fHistNP_shuffleV0M"; + AliTHn *fHistNP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistNP) { + Printf("fHistNP %s not found!!!",histoName.Data()); + break; + } + fHistNP->FillParent(); fHistNP->DeleteContainers(); + + if(!kShuffling) + histoName = "fHistPPV0M"; + else if(kShuffling) + histoName = "fHistPP_shuffleV0M"; + AliTHn *fHistPP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistPP) { + Printf("fHistPP %s not found!!!",histoName.Data()); + break; + } + fHistPP->FillParent(); fHistPP->DeleteContainers(); + + if(!kShuffling) + histoName = "fHistNNV0M"; + else if(kShuffling) + histoName = "fHistNN_shuffleV0M"; + AliTHn *fHistNN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistNN) { + Printf("fHistNN %s not found!!!",histoName.Data()); + break; + } + fHistNN->FillParent(); fHistNN->DeleteContainers(); + + return listBF; +} + +//______________________________________________________// +void draw(TList *listBF, TList *listBFShuffled, + Double_t psiMin, Double_t psiMax) { + gROOT->LoadMacro("~/SetPlotStyle.C"); + SetPlotStyle(); + gStyle->SetPalette(1,0); + + //balance function + AliTHn *hP = NULL; + AliTHn *hN = NULL; + AliTHn *hPN = NULL; + AliTHn *hNP = NULL; + AliTHn *hPP = NULL; + AliTHn *hNN = NULL; + //listBF->ls(); + //Printf("================="); + hP = (AliTHn*) listBF->FindObject("fHistPV0M"); + hN = (AliTHn*) listBF->FindObject("fHistNV0M"); + hPN = (AliTHn*) listBF->FindObject("fHistPNV0M"); + hNP = (AliTHn*) listBF->FindObject("fHistNPV0M"); + hPP = (AliTHn*) listBF->FindObject("fHistPPV0M"); + hNN = (AliTHn*) listBF->FindObject("fHistNNV0M"); + + AliBalancePsi *b = new AliBalancePsi(); + b->SetHistNp(hP); + b->SetHistNn(hN); + b->SetHistNpn(hPN); + b->SetHistNnp(hNP); + b->SetHistNpp(hPP); + b->SetHistNnn(hNN); + + //balance function shuffling + AliTHn *hPShuffled = NULL; + AliTHn *hNShuffled = NULL; + AliTHn *hPNShuffled = NULL; + AliTHn *hNPShuffled = NULL; + AliTHn *hPPShuffled = NULL; + AliTHn *hNNShuffled = NULL; + //listBFShuffled->ls(); + hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M"); + hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M"); + hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M"); + hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M"); + hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M"); + hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M"); + + AliBalancePsi *bShuffled = new AliBalancePsi(); + bShuffled->SetHistNp(hPShuffled); + bShuffled->SetHistNn(hNShuffled); + bShuffled->SetHistNpn(hPNShuffled); + bShuffled->SetHistNnp(hNPShuffled); + bShuffled->SetHistNpp(hPPShuffled); + bShuffled->SetHistNnn(hNNShuffled); + + TH1D *gHistBalanceFunction[numberOfCentralityBins]; + TH1D *gHistBalanceFunctionShuffled[numberOfCentralityBins]; + TCanvas *c1[numberOfCentralityBins]; + TString histoTitle, pngName; + TLegend *legend[numberOfCentralityBins]; + + //loop over the centrality bins + for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) { + histoTitle = "Centrality: "; + histoTitle += centralityArray[iCentralityBin]; + histoTitle += "%"; + histoTitle += " | "; histoTitle += psiMin; + histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax; + + gHistBalanceFunction[iCentralityBin] = b->GetBalanceFunctionHistogram(2,2,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + //gHistBalanceFunction[iCentralityBin] = b->GetBalanceFunctionHistogram(3,3,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistBalanceFunction[iCentralityBin]->SetMarkerStyle(20); + gHistBalanceFunction[iCentralityBin]->SetTitle(histoTitle.Data()); + gHistBalanceFunction[iCentralityBin]->GetYaxis()->SetTitleOffset(1.3); + + gHistBalanceFunctionShuffled[iCentralityBin] = bShuffled->GetBalanceFunctionHistogram(2,2,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + //gHistBalanceFunctionShuffled[iCentralityBin] = b->GetBalanceFunctionHistogram(3,3,gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistBalanceFunctionShuffled[iCentralityBin]->SetMarkerStyle(24); + + c1[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,0,600,500); + c1[iCentralityBin]->SetFillColor(10); + c1[iCentralityBin]->SetHighLightColor(10); + c1[iCentralityBin]->SetLeftMargin(0.15); + gHistBalanceFunction[iCentralityBin]->Draw("E"); + gHistBalanceFunctionShuffled[iCentralityBin]->Draw("ESAME"); + + legend[iCentralityBin] = new TLegend(0.18,0.6,0.45,0.82,"","brNDC"); + legend[iCentralityBin]->SetTextSize(0.045); + legend[iCentralityBin]->SetTextFont(42); + legend[iCentralityBin]->SetBorderSize(0); + legend[iCentralityBin]->SetFillStyle(0); + legend[iCentralityBin]->SetFillColor(10); + legend[iCentralityBin]->SetMargin(0.25); + legend[iCentralityBin]->SetShadowColor(10); + legend[iCentralityBin]->AddEntry(gHistBalanceFunction[iCentralityBin],"Data","lp"); + legend[iCentralityBin]->AddEntry(gHistBalanceFunctionShuffled[iCentralityBin],"Shuffled data","lp"); + legend[iCentralityBin]->Draw(); + + pngName = "BalanceFunctionDeltaEta.Centrality"; + pngName += centralityArray[iCentralityBin]; + pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax; + pngName += ".png"; + c1[iCentralityBin]->SaveAs(pngName.Data()); + + GetWeightedMean(gHistBalanceFunction[iCentralityBin]); + GetWeightedMean(gHistBalanceFunctionShuffled[iCentralityBin]); + + TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex; + meanLatex = "#mu = "; + meanLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetMean()); + meanLatex += " #pm "; + meanLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetMeanError()); + + rmsLatex = "#sigma = "; + rmsLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetRMS()); + rmsLatex += " #pm "; + rmsLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetRMSError()); + + skewnessLatex = "S = "; + skewnessLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetSkewness(1)); + skewnessLatex += " #pm "; + skewnessLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetSkewness(11)); + + kurtosisLatex = "K = "; + kurtosisLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetKurtosis(1)); + kurtosisLatex += " #pm "; + kurtosisLatex += Form("%.3f",gHistBalanceFunction[iCentralityBin]->GetKurtosis(11)); + + TLatex *latex = new TLatex(); + latex->SetNDC(); + latex->SetTextSize(0.035); + latex->SetTextColor(1); + latex->DrawLatex(0.64,0.85,meanLatex.Data()); + latex->DrawLatex(0.64,0.81,rmsLatex.Data()); + latex->DrawLatex(0.64,0.77,skewnessLatex.Data()); + latex->DrawLatex(0.64,0.73,kurtosisLatex.Data()); + Printf("Mean: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetMean(),gHistBalanceFunction[iCentralityBin]->GetMeanError()); + Printf("RMS: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetRMS(),gHistBalanceFunction[iCentralityBin]->GetRMSError()); + Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetSkewness(1),gHistBalanceFunction[iCentralityBin]->GetSkewness(11)); + Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunction[iCentralityBin]->GetKurtosis(1),gHistBalanceFunction[iCentralityBin]->GetKurtosis(11)); + } +} + +//____________________________________________________________________// +void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) { + //Prints the calculated width of the BF and its error + Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0; + Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0; + Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0; + Double_t deltaBalP2 = 0.0, integral = 0.0; + Double_t deltaErrorNew = 0.0; + + //Retrieve this variables from Histogram + Int_t fNumberOfBins = gHistBalance->GetNbinsX(); + Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning! + + cout<<"=================================================="<GetBinContent(i)<<"\t Error: "<GetBinError(i)<<"\t bin: "<GetBinCenter(i))<GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi| + gSumBi += gHistBalance->GetBinContent(i); + gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i)); + gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2); + gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2); + gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2); + gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2); + + deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2); + integral += fP2Step*gHistBalance->GetBinContent(i); + } + for(Int_t i = fStartBin; i < fNumberOfBins; i++) + deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2); + + Double_t integralError = TMath::Sqrt(deltaBalP2); + + Double_t delta = gSumBiXi / gSumBi; + Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) ); + cout<<"=================================================="<Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + gSystem->Load("libEventMixing.so"); + gSystem->Load("libCORRFW.so"); + gSystem->Load("libPWGTools.so"); + gSystem->Load("libPWGCFebye.so"); + + //Prepare the objects and return them + TList *list = GetListOfObjects(filename); + if(!list) { + Printf("The TList object was not created"); + return; + } + else + draw(list,psiMin,psiMax); +} + +//______________________________________________________// +TList *GetListOfObjects(const char* filename) { + //Get the TList objects (QA, bf, bf shuffled) + TList *listQA = 0x0; + TList *listBF = 0x0; + TList *listBFShuffling = 0x0; + + //Open the file + TFile *f = TFile::Open(filename); + if((!f)||(!f->IsOpen())) { + Printf("The file %s is not found. Aborting...",filename); + return listBF; + } + //f->ls(); + + TDirectoryFile *dir = dynamic_cast(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis")); + if(!dir) { + Printf("The TDirectoryFile is not found. Aborting...",filename); + return listBF; + } + //dir->ls(); + + TString listBFName = "listBF_0-100_V0M_vZ10.0_DCAxy-1.0_DCAz-1.0_Pt0.3-5.0_Eta-0.8-0.8_Chi-1.0_nClus-1_Bit1_withCentralTrigger"; + listBF = dynamic_cast(dir->Get(listBFName.Data())); + //listBF->ls(); + + //Get the histograms + TString histoName = "fHistPV0M"; + AliTHn *fHistP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistP) { + Printf("fHistP %s not found!!!",histoName.Data()); + break; + } + fHistP->FillParent(); fHistP->DeleteContainers(); + + histoName = "fHistNV0M"; + AliTHn *fHistN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistN) { + Printf("fHistN %s not found!!!",histoName.Data()); + break; + } + fHistN->FillParent(); fHistN->DeleteContainers(); + + histoName = "fHistPNV0M"; + AliTHn *fHistPN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistPN) { + Printf("fHistPN %s not found!!!",histoName.Data()); + break; + } + fHistPN->FillParent(); fHistPN->DeleteContainers(); + + histoName = "fHistNPV0M"; + AliTHn *fHistNP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistNP) { + Printf("fHistNP %s not found!!!",histoName.Data()); + break; + } + fHistNP->FillParent(); fHistNP->DeleteContainers(); + + histoName = "fHistPPV0M"; + AliTHn *fHistPP = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistPP) { + Printf("fHistPP %s not found!!!",histoName.Data()); + break; + } + fHistPP->FillParent(); fHistPP->DeleteContainers(); + + histoName = "fHistNNV0M"; + AliTHn *fHistNN = dynamic_cast(listBF->FindObject(histoName.Data())); + if(!fHistNN) { + Printf("fHistNN %s not found!!!",histoName.Data()); + break; + } + fHistNN->FillParent(); fHistNN->DeleteContainers(); + + return listBF; +} + +//______________________________________________________// +void draw(TList *list, Double_t psiMin, Double_t psiMax) { + //Draws the correlation functions for every centrality bin + //(+-), (-+), (++), (--) + gROOT->LoadMacro("~/SetPlotStyle.C"); + SetPlotStyle(); + gStyle->SetPalette(1,0); + + AliTHn *hP = NULL; + AliTHn *hN = NULL; + AliTHn *hPN = NULL; + AliTHn *hNP = NULL; + AliTHn *hPP = NULL; + AliTHn *hNN = NULL; + + hP = (AliTHn*) list->FindObject("fHistPV0M"); + hN = (AliTHn*) list->FindObject("fHistNV0M"); + hPN = (AliTHn*) list->FindObject("fHistPNV0M"); + hNP = (AliTHn*) list->FindObject("fHistNPV0M"); + hPP = (AliTHn*) list->FindObject("fHistPPV0M"); + hNN = (AliTHn*) list->FindObject("fHistNNV0M"); + + //Create the AliBalancePsi object and fill it with the AliTHn objects + AliBalancePsi *b = new AliBalancePsi(); + b->SetHistNp(hP); + b->SetHistNn(hN); + b->SetHistNpn(hPN); + b->SetHistNnp(hNP); + b->SetHistNpp(hPP); + b->SetHistNnn(hNN); + TH2D *gHistPN[numberOfCentralityBins]; + TH2D *gHistNP[numberOfCentralityBins]; + TH2D *gHistPP[numberOfCentralityBins]; + TH2D *gHistNN[numberOfCentralityBins]; + + TCanvas *cPN[numberOfCentralityBins]; + TCanvas *cNP[numberOfCentralityBins]; + TCanvas *cPP[numberOfCentralityBins]; + TCanvas *cNN[numberOfCentralityBins]; + TString histoTitle, pngName; + //loop over the centrality bins + for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) { + histoTitle = "(+-) | Centrality: "; + histoTitle += centralityArray[iCentralityBin]; + histoTitle += "%"; + histoTitle += " | "; histoTitle += psiMin; + histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax; + gHistPN[iCentralityBin] = b->GetCorrelationFunctionPN(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistPN[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5); + gHistPN[iCentralityBin]->SetTitle(histoTitle.Data()); + cPN[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,0,400,400); + cPN[iCentralityBin]->SetFillColor(10); + cPN[iCentralityBin]->SetHighLightColor(10); + gHistPN[iCentralityBin]->DrawCopy("lego2"); + pngName = "DeltaPhiDeltaEta.Centrality"; + pngName += centralityArray[iCentralityBin]; + pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax; + pngName += ".PositiveNegative.png"; + cPN[iCentralityBin]->SaveAs(pngName.Data()); + + histoTitle = "(-+) | Centrality: "; + histoTitle += centralityArray[iCentralityBin]; + histoTitle += "%"; + histoTitle += " | "; histoTitle += psiMin; + histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax; + gHistNP[iCentralityBin] = b->GetCorrelationFunctionNP(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistNP[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5); + gHistNP[iCentralityBin]->SetTitle(histoTitle.Data()); + cNP[iCentralityBin] = new TCanvas(histoTitle.Data(),"",400,0,400,400); + cNP[iCentralityBin]->SetFillColor(10); + cNP[iCentralityBin]->SetHighLightColor(10); + gHistNP[iCentralityBin]->DrawCopy("lego2"); + pngName = "DeltaPhiDeltaEta.Centrality"; + pngName += centralityArray[iCentralityBin]; + pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax; + pngName += ".NegativePositive.png"; + cNP[iCentralityBin]->SaveAs(pngName.Data()); + + histoTitle = "(++) | Centrality: "; + histoTitle += centralityArray[iCentralityBin]; + histoTitle += "%"; + histoTitle += " | "; histoTitle += psiMin; + histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax; + gHistPP[iCentralityBin] = b->GetCorrelationFunctionPP(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistPP[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5); + gHistPP[iCentralityBin]->SetTitle(histoTitle.Data()); + cPP[iCentralityBin] = new TCanvas(histoTitle.Data(),"",0,400,400,400); + cPP[iCentralityBin]->SetFillColor(10); + cPP[iCentralityBin]->SetHighLightColor(10); + gHistPP[iCentralityBin]->DrawCopy("lego2"); + pngName = "DeltaPhiDeltaEta.Centrality"; + pngName += centralityArray[iCentralityBin]; + pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax; + pngName += ".PositivePositive.png"; + cPP[iCentralityBin]->SaveAs(pngName.Data()); + + histoTitle = "(--) | Centrality: "; + histoTitle += centralityArray[iCentralityBin]; + histoTitle += "%"; + histoTitle += " | "; histoTitle += psiMin; + histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax; + gHistNN[iCentralityBin] = b->GetCorrelationFunctionNN(gMinCentrality[iCentralityBin],gMaxCentrality[iCentralityBin],psiMin,psiMax); + gHistNN[iCentralityBin]->GetYaxis()->SetTitleOffset(1.5); + gHistNN[iCentralityBin]->SetTitle(histoTitle.Data()); + cNN[iCentralityBin] = new TCanvas(histoTitle.Data(),"",400,400,400,400); + cNN[iCentralityBin]->SetFillColor(10); + cNN[iCentralityBin]->SetHighLightColor(10); + gHistNN[iCentralityBin]->DrawCopy("lego2"); + pngName = "DeltaPhiDeltaEta.Centrality"; + pngName += centralityArray[iCentralityBin]; + pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax; + pngName += ".NegativeNegative.png"; + cNN[iCentralityBin]->SaveAs(pngName.Data()); + }//end of loop over centralities +} + -- 2.43.0