-const Int_t numberOfCentralityBins = 8;
-TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
+const Int_t numberOfCentralityBins = 12;
+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;
void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
Int_t gCentrality = 1,
+ Int_t gBit = -1,
+ const char* gCentralityEstimator = 0x0,
+ Bool_t kShowShuffled = kFALSE,
+ Bool_t kShowMixed = kTRUE,
Double_t psiMin = -0.5, Double_t psiMax = 0.5,
+ Double_t vertexZMin = -10.,
+ Double_t vertexZMax = 10.,
Double_t ptTriggerMin = -1.,
Double_t ptTriggerMax = -1.,
Double_t ptAssociatedMin = -1.,
- Double_t ptAssociatedMax = -1.) {
+ Double_t ptAssociatedMax = -1.,
+ Bool_t kUseVzBinning = kTRUE,
+ Bool_t k2pMethod = kTRUE,
+ TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
+ Bool_t bToy = kFALSE)
+{
//Macro that draws the BF distributions for each centrality bin
//for reaction plane dependent analysis
//Author: Panos.Christakoglou@nikhef.nl
gSystem->Load("libPWGTools.so");
gSystem->Load("libPWGCFebye.so");
+ //gROOT->LoadMacro("~/SetPlotStyle.C");
+ //SetPlotStyle();
+ gStyle->SetPalette(1,0);
+
//Prepare the objects and return them
- TList *listBF = GetListOfObjects(filename,gCentrality,0);
- TList *listBFShuffled = GetListOfObjects(filename,gCentrality,1);
- TList *listBFMixed = GetListOfObjects(filename,gCentrality,2);
+ TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy);
+ TList *listBFShuffled = NULL;
+ if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy);
+ TList *listBFMixed = NULL;
+ if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy);
if(!listBF) {
Printf("The TList object was not created");
return;
}
else
- draw(listBF,listBFShuffled,listBFMixed,gCentrality,psiMin,psiMax,
- ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+ draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
+ psiMin,psiMax,vertexZMin,vertexZMax,
+ ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
+ kUseVzBinning,k2pMethod,eventClass);
}
//______________________________________________________//
-TList *GetListOfObjects(const char* filename,
- Int_t gCentrality, Int_t kData = 1) {
+TList *GetListOfObjects(const char* filename,
+ Int_t gCentrality,
+ Int_t gBit,
+ const char *gCentralityEstimator,
+ Int_t kData = 1,
+ Bool_t bToy = kFALSE) {
//Get the TList objects (QA, bf, bf shuffled)
- //kData == 0: data
- //kData == 1: shuffling
- //kData == 2: mixing
- TList *listQA = 0x0;
TList *listBF = 0x0;
//Open the file
- TFile *f = TFile::Open(filename);
+ TFile *f = TFile::Open(filename,"UPDATE");
if((!f)||(!f->IsOpen())) {
Printf("The file %s is not found. Aborting...",filename);
return listBF;
TString listBFName;
if(kData == 0) {
//cout<<"no shuffling - no mixing"<<endl;
- listBFName = "listBFPsi_";
+ listBFName = "listBFPsi";
}
else if(kData == 1) {
//cout<<"shuffling - no mixing"<<endl;
- listBFName = "listBFPsiShuffled_";
+ listBFName = "listBFPsiShuffled";
}
else if(kData == 2) {
//cout<<"no shuffling - mixing"<<endl;
- listBFName = "listBFPsiMixed_";
+ listBFName = "listBFPsiMixed";
}
- listBFName += centralityArray[gCentrality-1];
- listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
- cout<<"======================================================="<<endl;
- cout<<"List name: "<<listBF->GetName()<<endl;
- //listBF->ls();
- //Get the histograms
- TString histoName;
- if(kData == 0)
- histoName = "fHistPV0M";
- else if(kData == 1)
- histoName = "fHistP_shuffleV0M";
- else if(kData == 2)
- histoName = "fHistPV0M";
- AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistP) {
- Printf("fHistP %s not found!!!",histoName.Data());
- break;
- }
- fHistP->FillParent(); fHistP->DeleteContainers();
-
- if(kData == 0)
- histoName = "fHistNV0M";
- if(kData == 1)
- histoName = "fHistN_shuffleV0M";
- if(kData == 2)
- histoName = "fHistNV0M";
- AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistN) {
- Printf("fHistN %s not found!!!",histoName.Data());
- break;
- }
- fHistN->FillParent(); fHistN->DeleteContainers();
+ // different list names in case of toy model
+ if(!bToy){
+ listBFName += "_"
+ listBFName += centralityArray[gCentrality-1];
+ if(gBit > -1) {
+ listBFName += "_Bit"; listBFName += gBit; }
+ if(gCentralityEstimator) {
+ listBFName += "_"; listBFName += gCentralityEstimator;}
+ }
+ else{
+ listBFName.ReplaceAll("Psi","");
+ }
+
+ // histograms were already retrieved (in first iteration)
+ if(dir->Get(Form("%s_histograms",listBFName.Data()))){
+ //cout<<"second iteration"<<endl;
+ listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
+ }
+
+ // histograms were not yet retrieved (this is the first iteration)
+ else{
+
+ listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
+ cout<<"======================================================="<<endl;
+ cout<<"List name: "<<listBF->GetName()<<endl;
+ //listBF->ls();
+
+ //Get the histograms
+ TString histoName;
+ if(kData == 0)
+ histoName = "fHistP";
+ else if(kData == 1)
+ histoName = "fHistP_shuffle";
+ else if(kData == 2)
+ histoName = "fHistP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistP) {
+ Printf("fHistP %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistP->FillParent(); fHistP->DeleteContainers();
+
+ if(kData == 0)
+ histoName = "fHistN";
+ if(kData == 1)
+ histoName = "fHistN_shuffle";
+ if(kData == 2)
+ histoName = "fHistN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistN) {
+ Printf("fHistN %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistN->FillParent(); fHistN->DeleteContainers();
+
+ if(kData == 0)
+ histoName = "fHistPN";
+ if(kData == 1)
+ histoName = "fHistPN_shuffle";
+ if(kData == 2)
+ histoName = "fHistPN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistPN) {
+ Printf("fHistPN %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistPN->FillParent(); fHistPN->DeleteContainers();
+
+ if(kData == 0)
+ histoName = "fHistNP";
+ if(kData == 1)
+ histoName = "fHistNP_shuffle";
+ if(kData == 2)
+ histoName = "fHistNP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistNP) {
+ Printf("fHistNP %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistNP->FillParent(); fHistNP->DeleteContainers();
+
+ if(kData == 0)
+ histoName = "fHistPP";
+ if(kData == 1)
+ histoName = "fHistPP_shuffle";
+ if(kData == 2)
+ histoName = "fHistPP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistPP) {
+ Printf("fHistPP %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistPP->FillParent(); fHistPP->DeleteContainers();
+
+ if(kData == 0)
+ histoName = "fHistNN";
+ if(kData == 1)
+ histoName = "fHistNN_shuffle";
+ if(kData == 2)
+ histoName = "fHistNN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
+ if(!fHistNN) {
+ Printf("fHistNN %s not found!!!",histoName.Data());
+ break;
+ }
+ fHistNN->FillParent(); fHistNN->DeleteContainers();
- if(kData == 0)
- histoName = "fHistPNV0M";
- if(kData == 1)
- histoName = "fHistPN_shuffleV0M";
- if(kData == 2)
- histoName = "fHistPNV0M";
- AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistPN) {
- Printf("fHistPN %s not found!!!",histoName.Data());
- break;
- }
- fHistPN->FillParent(); fHistPN->DeleteContainers();
+ dir->cd();
+ listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
+
+ }// first iteration
- if(kData == 0)
- histoName = "fHistNPV0M";
- if(kData == 1)
- histoName = "fHistNP_shuffleV0M";
- if(kData == 2)
- histoName = "fHistNPV0M";
- AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistNP) {
- Printf("fHistNP %s not found!!!",histoName.Data());
- break;
- }
- fHistNP->FillParent(); fHistNP->DeleteContainers();
-
- if(kData == 0)
- histoName = "fHistPPV0M";
- if(kData == 1)
- histoName = "fHistPP_shuffleV0M";
- if(kData == 2)
- histoName = "fHistPPV0M";
- AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistPP) {
- Printf("fHistPP %s not found!!!",histoName.Data());
- break;
- }
- fHistPP->FillParent(); fHistPP->DeleteContainers();
-
- if(kData == 0)
- histoName = "fHistNNV0M";
- if(kData == 1)
- histoName = "fHistNN_shuffleV0M";
- if(kData == 2)
- histoName = "fHistNNV0M";
- AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
- if(!fHistNN) {
- Printf("fHistNN %s not found!!!",histoName.Data());
- break;
- }
- fHistNN->FillParent(); fHistNN->DeleteContainers();
+ f->Close();
return listBF;
}
//______________________________________________________//
void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
- Int_t gCentrality, Double_t psiMin, Double_t psiMax,
+ Int_t gCentrality, const char* gCentralityEstimator,
+ Double_t psiMin, Double_t psiMax,
+ Double_t vertexZMin,
+ Double_t vertexZMax,
Double_t ptTriggerMin, Double_t ptTriggerMax,
- Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
- gROOT->LoadMacro("~/SetPlotStyle.C");
- SetPlotStyle();
- gStyle->SetPalette(1,0);
-
+ Double_t ptAssociatedMin, Double_t ptAssociatedMax,
+ Bool_t kUseVzBinning=kFALSE,
+ Bool_t k2pMethod = kFALSE, TString eventClass) {
//balance function
AliTHn *hP = NULL;
AliTHn *hN = 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");
+ TString histoName = "fHistP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hP = (AliTHn*) listBF->FindObject(histoName.Data());
+ hP->SetName("gHistP");
+ histoName = "fHistN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hN = (AliTHn*) listBF->FindObject(histoName.Data());
+ hN->SetName("gHistN");
+ histoName = "fHistPN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPN = (AliTHn*) listBF->FindObject(histoName.Data());
+ hPN->SetName("gHistPN");
+ histoName = "fHistNP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNP = (AliTHn*) listBF->FindObject(histoName.Data());
+ hNP->SetName("gHistNP");
+ histoName = "fHistPP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPP = (AliTHn*) listBF->FindObject(histoName.Data());
+ hPP->SetName("gHistPP");
+ histoName = "fHistNN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNN = (AliTHn*) listBF->FindObject(histoName.Data());
+ hNN->SetName("gHistNN");
AliBalancePsi *b = new AliBalancePsi();
+ b->SetEventClass(eventClass);
b->SetHistNp(hP);
b->SetHistNn(hN);
b->SetHistNpn(hPN);
b->SetHistNnp(hNP);
b->SetHistNpp(hPP);
b->SetHistNnn(hNN);
+ if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
+
//balance function shuffling
AliTHn *hPShuffled = 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);
+ if(listBFShuffled) {
+ //listBFShuffled->ls();
+ histoName = "fHistP_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hPShuffled->SetName("gHistPShuffled");
+ histoName = "fHistN_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hNShuffled->SetName("gHistNShuffled");
+ histoName = "fHistPN_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hPNShuffled->SetName("gHistPNShuffled");
+ histoName = "fHistNP_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hNPShuffled->SetName("gHistNPShuffled");
+ histoName = "fHistPP_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hPPShuffled->SetName("gHistPPShuffled");
+ histoName = "fHistNN_shuffle";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
+ hNNShuffled->SetName("gHistNNShuffled");
+
+ AliBalancePsi *bShuffled = new AliBalancePsi();
+ bShuffled->SetEventClass(eventClass);
+ bShuffled->SetHistNp(hPShuffled);
+ bShuffled->SetHistNn(hNShuffled);
+ bShuffled->SetHistNpn(hPNShuffled);
+ bShuffled->SetHistNnp(hNPShuffled);
+ bShuffled->SetHistNpp(hPPShuffled);
+ bShuffled->SetHistNnn(hNNShuffled);
+ if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
+
+ }
//balance function mixing
AliTHn *hPMixed = NULL;
AliTHn *hNPMixed = NULL;
AliTHn *hPPMixed = NULL;
AliTHn *hNNMixed = NULL;
- //listBFMixed->ls();
- hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
- hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
- hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
- hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
- hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
- hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
-
- AliBalancePsi *bMixed = new AliBalancePsi();
- bMixed->SetHistNp(hPMixed);
- bMixed->SetHistNn(hNMixed);
- bMixed->SetHistNpn(hPNMixed);
- bMixed->SetHistNnp(hNPMixed);
- bMixed->SetHistNpp(hPPMixed);
- bMixed->SetHistNnn(hNNMixed);
+
+ if(listBFMixed) {
+ //listBFMixed->ls();
+ histoName = "fHistP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ hPMixed->SetName("gHistPMixed");
+ histoName = "fHistN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ hNMixed->SetName("gHistNMixed");
+ histoName = "fHistPN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ hPNMixed->SetName("gHistPNMixed");
+ histoName = "fHistNP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ histoName = "fHistNP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNPMixed->SetName("gHistNPMixed");
+ histoName = "fHistPP";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ hPPMixed->SetName("gHistPPMixed");
+ histoName = "fHistNN";
+ if(gCentralityEstimator) histoName += gCentralityEstimator;
+ hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
+ hNNMixed->SetName("gHistNNMixed");
+
+ AliBalancePsi *bMixed = new AliBalancePsi();
+ bMixed->SetEventClass(eventClass);
+ bMixed->SetHistNp(hPMixed);
+ bMixed->SetHistNn(hNMixed);
+ bMixed->SetHistNpn(hPNMixed);
+ bMixed->SetHistNnp(hNPMixed);
+ bMixed->SetHistNpp(hPPMixed);
+ bMixed->SetHistNnn(hNNMixed);
+ if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
+
+ }
TH2D *gHistBalanceFunction;
TH2D *gHistBalanceFunctionSubtracted;
TH2D *gHistBalanceFunctionMixed;
TString histoTitle, pngName;
- histoTitle = "Centrality: ";
- histoTitle += centralityArray[gCentrality-1];
- histoTitle += "%";
- if((psiMin == -0.5)&&(psiMax == 0.5))
- histoTitle += " (-7.5^{o} < #phi - #Psi_{2} < 7.5^{o})";
- else if((psiMin == 0.5)&&(psiMax == 1.5))
- histoTitle += " (37.5^{o} < #phi - #Psi_{2} < 52.5^{o})";
- else if((psiMin == 1.5)&&(psiMax == 2.5))
- histoTitle += " (82.5^{o} < #phi - #Psi_{2} < 97.5^{o})";
- else
- histoTitle += " (0^{o} < #phi - #Psi_{2} < 180^{o})";
-
- gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+ if(eventClass == "Centrality"){
+ histoTitle = "Centrality: ";
+ histoTitle += psiMin;
+ histoTitle += " - ";
+ histoTitle += psiMax;
+ histoTitle += " % ";
+ histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
+ }
+ else if(eventClass == "Multiplicity"){
+ histoTitle = "Multiplicity: ";
+ histoTitle += psiMin;
+ histoTitle += " - ";
+ histoTitle += psiMax;
+ histoTitle += " tracks";
+ histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
+ }
+ else{ // "EventPlane" (default)
+ histoTitle = "Centrality: ";
+ histoTitle += centralityArray[gCentrality-1];
+ histoTitle += "%";
+ if((psiMin == -0.5)&&(psiMax == 0.5))
+ histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
+ else if((psiMin == 0.5)&&(psiMax == 1.5))
+ histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
+ else if((psiMin == 1.5)&&(psiMax == 2.5))
+ histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
+ else
+ histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
+ }
+
+ if(k2pMethod)
+ if(bMixed)
+ gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
+ else{
+ cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
+ return;
+ }
+ else
+ gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
gHistBalanceFunction->SetTitle(histoTitle.Data());
gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
gHistBalanceFunction->SetName("gHistBalanceFunction");
- gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
- gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
- gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
-
- gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
- gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
- gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
+ if(listBFShuffled) {
+
+ if(k2pMethod)
+ if(bMixed)
+ gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
+ else{
+ cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
+ return;
+ }
+ else
+ gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+ gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
+ gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
+ }
- gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
- gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
- gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
- gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
+ if(listBFMixed) {
+ if(k2pMethod)
+ if(bMixed)
+ gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
+ else{
+ cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
+ return;
+ }
+ else
+ gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+ gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
+ gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
+
+ gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
+ gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
+ gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
+ }
//Draw the results
TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
c1->SetHighLightColor(10);
c1->SetLeftMargin(0.15);
gHistBalanceFunction->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
c1a->SetFillColor(10);
c1a->SetHighLightColor(10);
c1a->SetLeftMargin(0.15);
gHistBalanceFunction->DrawCopy("colz");
- TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
- c2->SetFillColor(10);
- c2->SetHighLightColor(10);
- c2->SetLeftMargin(0.15);
- gHistBalanceFunctionShuffled->DrawCopy("lego2");
- TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
- c2a->SetFillColor(10);
- c2a->SetHighLightColor(10);
- c2a->SetLeftMargin(0.15);
- gHistBalanceFunctionShuffled->DrawCopy("colz");
-
- TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
- c3->SetFillColor(10);
- c3->SetHighLightColor(10);
- c3->SetLeftMargin(0.15);
- gHistBalanceFunctionMixed->DrawCopy("lego2");
- TCanvas *c3a = new TCanvas("c3a","",700,200,600,500);
- c3a->SetFillColor(10);
- c3a->SetHighLightColor(10);
- c3a->SetLeftMargin(0.15);
- gHistBalanceFunctionMixed->DrawCopy("colz");
-
- TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
- c4->SetFillColor(10);
- c4->SetHighLightColor(10);
- c4->SetLeftMargin(0.15);
- gHistBalanceFunctionSubtracted->DrawCopy("lego2");
- TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
- c4a->SetFillColor(10);
- c4a->SetHighLightColor(10);
- c4a->SetLeftMargin(0.15);
- gHistBalanceFunctionSubtracted->DrawCopy("colz");
-
- TString newFileName = "balanceFunction2D.Centrality";
- newFileName += gCentrality; newFileName += ".Psi";
- if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
- else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
- else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
- else newFileName += "0.PttFrom";
- newFileName += ptTriggerMin; newFileName += "To";
- newFileName += ptTriggerMax; newFileName += ".PtaFrom";
- newFileName += ptAssociatedMin; newFileName += "To";
- newFileName += ptAssociatedMax; newFileName += ".root";
+ if(listBFShuffled) {
+ TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
+ c2->SetFillColor(10);
+ c2->SetHighLightColor(10);
+ c2->SetLeftMargin(0.15);
+ gHistBalanceFunctionShuffled->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+ TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
+ c2a->SetFillColor(10);
+ c2a->SetHighLightColor(10);
+ c2a->SetLeftMargin(0.15);
+ gHistBalanceFunctionShuffled->DrawCopy("colz");
+ }
+
+ if(listBFMixed) {
+ TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
+ c3->SetFillColor(10);
+ c3->SetHighLightColor(10);
+ c3->SetLeftMargin(0.15);
+ gHistBalanceFunctionMixed->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+ TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
+ c3a->SetFillColor(10);
+ c3a->SetHighLightColor(10);
+ c3a->SetLeftMargin(0.15);
+ gHistBalanceFunctionMixed->DrawCopy("colz");
+
+ TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
+ c4->SetFillColor(10);
+ c4->SetHighLightColor(10);
+ c4->SetLeftMargin(0.15);
+ gHistBalanceFunctionSubtracted->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+ TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
+ c4a->SetFillColor(10);
+ c4a->SetHighLightColor(10);
+ c4a->SetLeftMargin(0.15);
+ gHistBalanceFunctionSubtracted->DrawCopy("colz");
+
+ fitbalanceFunction(gCentrality, psiMin , psiMax,
+ ptTriggerMin, ptTriggerMax,
+ ptAssociatedMin, ptAssociatedMax,
+ gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
+ }
+
+ TString newFileName = "balanceFunction2D.";
+ if(eventClass == "Centrality"){
+ newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
+ newFileName += ".PsiAll.PttFrom";
+ }
+ else if(eventClass == "Multiplicity"){
+ newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
+ newFileName += ".PsiAll.PttFrom";
+ }
+ else{ // "EventPlane" (default)
+ newFileName += "Centrality";
+ newFileName += gCentrality; newFileName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
+ else newFileName += "All.PttFrom";
+ }
+ newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
+ newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
+ 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 *fOutput = new TFile(newFileName.Data(),"recreate");
fOutput->cd();
+ /*hP->Write(); hN->Write();
+ hPN->Write(); hNP->Write();
+ hPP->Write(); hNN->Write();
+ hPShuffled->Write(); hNShuffled->Write();
+ hPNShuffled->Write(); hNPShuffled->Write();
+ hPPShuffled->Write(); hNNShuffled->Write();
+ hPMixed->Write(); hNMixed->Write();
+ hPNMixed->Write(); hNPMixed->Write();
+ hPPMixed->Write(); hNNMixed->Write();*/
gHistBalanceFunction->Write();
- gHistBalanceFunctionShuffled->Write();
- gHistBalanceFunctionMixed->Write();
- gHistBalanceFunctionSubtracted->Write();
+ if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
+ if(listBFMixed) {
+ gHistBalanceFunctionMixed->Write();
+ gHistBalanceFunctionSubtracted->Write();
+ }
fOutput->Close();
}
//____________________________________________________________//
-void drawBFPsi2D(Int_t gCentrality = 1,
+void fitbalanceFunction(Int_t gCentrality = 1,
+ Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+ Double_t ptTriggerMin = -1.,
+ Double_t ptTriggerMax = -1.,
+ Double_t ptAssociatedMin = -1.,
+ Double_t ptAssociatedMax = -1.,
+ TH2D *gHist,
+ Bool_t k2pMethod = kFALSE,
+ TString eventClass="EventPlane") {
+ //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
+ //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
+ cout<<"FITTING FUNCTION"<<endl;
+
+ TF2 *gFitFunction = new TF2("gFitFunction","[0] + [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2)) + [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))",-1.2,1.2,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
+ gFitFunction->SetName("gFitFunction");
+
+ //Normalization
+ gFitFunction->SetParName(0,"N1");
+ gFitFunction->SetParameter(0,1.0);
+
+ //2D balance function
+ gFitFunction->SetParName(1,"N_{BF}");
+ gFitFunction->SetParameter(1,1.0);
+ gFitFunction->SetParLimits(1, 0., 100.);
+ gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
+ gFitFunction->SetParameter(2,0.6);
+ gFitFunction->SetParLimits(2, 0., 1.);
+ gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
+ gFitFunction->SetParameter(3,0.0);
+ gFitFunction->SetParLimits(3, -0.2, 0.2);
+ gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
+ gFitFunction->SetParameter(4,0.6);
+ gFitFunction->SetParLimits(4, 0., 1.);
+ gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
+ gFitFunction->SetParameter(5,0.0);
+ gFitFunction->SetParLimits(5, -0.2, 0.2);
+
+ //short range structure
+ gFitFunction->SetParName(6,"N_{SR}");
+ gFitFunction->SetParameter(6,5.0);
+ gFitFunction->SetParLimits(6, 0., 100.);
+ gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
+ gFitFunction->SetParameter(7,0.01);
+ gFitFunction->SetParLimits(7, 0.0, 0.1);
+ gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
+ gFitFunction->SetParameter(8,0.0);
+ gFitFunction->SetParLimits(8, -0.01, 0.01);
+ gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
+ gFitFunction->SetParameter(9,0.01);
+ gFitFunction->SetParLimits(9, 0.0, 0.1);
+ gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
+ gFitFunction->SetParameter(10,0.0);
+ gFitFunction->SetParLimits(10, -0.01, 0.01);
+
+
+ //Cloning the histogram
+ TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
+ gHistResidual->SetName("gHistResidual");
+ gHistResidual->Sumw2();
+
+ //Fitting the 2D bf
+ for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
+ gHist->Fit("gFitFunction","nm");
+ for(Int_t iParam = 0; iParam < 11; iParam++)
+ gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
+ }
+ cout<<"======================================================"<<endl;
+ cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
+ cout<<"======================================================"<<endl;
+
+ //Getting the residual
+ gHistResidual->Add(gFitFunction,-1);
+
+ //Write to output file
+ TString newFileName = "balanceFunctionFit2D.";
+ if(eventClass == "Centrality"){
+ newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
+ newFileName += ".PsiAll.PttFrom";
+ }
+ else if(eventClass == "Multiplicity"){
+ newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
+ newFileName += ".PsiAll.PttFrom";
+ }
+ else{ // "EventPlane" (default)
+ newFileName += "Centrality";
+ newFileName += gCentrality; newFileName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
+ else newFileName += "All.PttFrom";
+ }
+ newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
+ newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
+ 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();
+ gFitFunction->Write();
+ newFile->Close();
+}
+
+//____________________________________________________________//
+void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
+ const char* gCentralityEstimator = "V0M",
+ Int_t gBit = 128,
+ const char* gEventPlaneEstimator = "VZERO",
+ Int_t gCentrality = 1,
+ Bool_t kShowShuffled = kFALSE,
+ Bool_t kShowMixed = kFALSE,
Double_t psiMin = -0.5, Double_t psiMax = 0.5,
Double_t ptTriggerMin = -1.,
Double_t ptTriggerMax = -1.,
Double_t ptAssociatedMin = -1.,
- Double_t ptAssociatedMax = -1.) {
+ Double_t ptAssociatedMax = -1.,
+ Bool_t k2pMethod = kTRUE) {
//Macro that draws the BF distributions for each centrality bin
//for reaction plane dependent analysis
//Author: Panos.Christakoglou@nikhef.nl
- gROOT->LoadMacro("~/SetPlotStyle.C");
- SetPlotStyle();
+ TGaxis::SetMaxDigits(3);
//Get the input file
- TString filename = "LHC11h/PttFrom";
- filename += ptTriggerMin; filename += "To";
- filename += ptTriggerMax; filename += "PtaFrom";
- filename += ptAssociatedMin; filename += "To";
- filename += ptAssociatedMax; filename += "/balanceFunction2D.Centrality";
+ TString filename = lhcPeriod;
+ filename += "/Centrality"; filename += gCentralityEstimator;
+ filename += "_Bit"; filename += gBit;
+ filename += "_"; filename += gEventPlaneEstimator;
+ filename +="/PttFrom";
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ filename += Form("%.1f",ptAssociatedMin); filename += "To";
+ filename += Form("%.1f",ptAssociatedMax);
+ filename += "/balanceFunction2D.Centrality";
filename += gCentrality; filename += ".Psi";
if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
- else filename += "0.PttFrom";
- filename += ptTriggerMin; filename += "To";
- filename += ptTriggerMax; filename += ".PtaFrom";
- filename += ptAssociatedMin; filename += "To";
- filename += ptAssociatedMax; filename += ".root";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
+ else filename += "All.PttFrom";
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ 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
TFile *f = TFile::Open(filename.Data());
gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
- gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
- gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #varphi)");
+ gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
+ gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
//Shuffled balance function
- TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
- gHistBalanceFunctionShuffled->SetStats(kFALSE);
- gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
- gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
- gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
- gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #varphi)");
+ if(kShowShuffled) {
+ TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
+ gHistBalanceFunctionShuffled->SetStats(kFALSE);
+ gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
+ gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
+ gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
+ gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
+ }
//Mixed balance function
- TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
- gHistBalanceFunctionMixed->SetStats(kFALSE);
- gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
- gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
- gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
- gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #varphi)");
+ if(kShowMixed) {
+ TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
+ gHistBalanceFunctionMixed->SetStats(kFALSE);
+ gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
+ gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
+ gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
+ gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
+ }
//Subtracted balance function
- TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
- gHistBalanceFunctionSubtracted->SetStats(kFALSE);
- gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
- gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
- gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
- gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
- gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #varphi)");
-
+ if(kShowMixed) {
+ TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
+ gHistBalanceFunctionSubtracted->SetStats(kFALSE);
+ gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
+ gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
+ gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
+ gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
+ }
+
TString pngName;
TString centralityLatex = "Centrality: ";
TLatex *latexInfo1 = new TLatex();
latexInfo1->SetNDC();
- latexInfo1->SetTextSize(0.04);
+ latexInfo1->SetTextSize(0.045);
latexInfo1->SetTextColor(1);
//Draw the results
TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
c1->SetFillColor(10); c1->SetHighLightColor(10);
c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
- gHistBalanceFunction->SetTitle("Raw balance function");
+ gHistBalanceFunction->SetTitle("");
gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
gHistBalanceFunction->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
+
+ TString pngName = "BalanceFunction2D.";
+ pngName += "Centrality";
+ pngName += gCentrality;
+ if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
+ else pngName += "All.PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ if(k2pMethod) pngName += "_2pMethod";
+ pngName += ".png";
+
+ c1->SaveAs(pngName.Data());
+
+ if(kShowShuffled) {
+ TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
+ c2->SetFillColor(10); c2->SetHighLightColor(10);
+ c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
+ gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
+ gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
+ gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionShuffled->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
- latexInfo1->DrawLatex(0.68,0.88,centralityLatex.Data());
- latexInfo1->DrawLatex(0.68,0.82,psiLatex.Data());
- latexInfo1->DrawLatex(0.68,0.76,pttLatex.Data());
- latexInfo1->DrawLatex(0.68,0.70,ptaLatex.Data());
-
- TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
- c2->SetFillColor(10); c2->SetHighLightColor(10);
- c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
- gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
- gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
- gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionShuffled->DrawCopy("lego2");
-
- latexInfo1->DrawLatex(0.68,0.88,centralityLatex.Data());
- latexInfo1->DrawLatex(0.68,0.82,psiLatex.Data());
- latexInfo1->DrawLatex(0.68,0.76,pttLatex.Data());
- latexInfo1->DrawLatex(0.68,0.70,ptaLatex.Data());
-
- TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
- c3->SetFillColor(10); c3->SetHighLightColor(10);
- c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
- gHistBalanceFunctionMixed->SetTitle("Mixed events");
- gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
- gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionMixed->DrawCopy("lego2");
+ latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
+ }
+
+ if(kShowMixed) {
+ TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
+ c3->SetFillColor(10); c3->SetHighLightColor(10);
+ c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
+ gHistBalanceFunctionMixed->SetTitle("Mixed events");
+ gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
+ gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionMixed->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
+ }
+
+ if(kShowMixed) {
+ TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
+ c4->SetFillColor(10); c4->SetHighLightColor(10);
+ c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
+ gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunctionSubtracted->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
+ latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
+ }
+}
+
+//____________________________________________________________//
+void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
+ Bool_t kProjectInEta = kFALSE,
+ Int_t binMin = 1,
+ Int_t binMax = 80,
+ Int_t gCentrality = 1,
+ Double_t psiMin = -0.5,
+ Double_t psiMax = 3.5,
+ Double_t ptTriggerMin = -1.,
+ Double_t ptTriggerMax = -1.,
+ Double_t ptAssociatedMin = -1.,
+ Double_t ptAssociatedMax = -1.,
+ Bool_t k2pMethod = kTRUE,
+ TString eventClass = "Centrality",
+ Bool_t bRootMoments = kFALSE) {
+ //Macro that draws the balance functions PROJECTIONS
+ //for each centrality bin for the different pT of trigger and
+ //associated particles
+ 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");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libPWGTools.so");
+ gSystem->Load("libPWGCFebye.so");
+
+ gStyle->SetOptStat(0);
+
+ //Get the input file
+ TString filename = "balanceFunction2D.";
+ if(eventClass == "Centrality"){
+ filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
+ filename += ".PsiAll.PttFrom";
+ }
+ else if(eventClass == "Multiplicity"){
+ filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
+ filename += ".PsiAll.PttFrom";
+ }
+ else{ // "EventPlane" (default)
+ filename += "Centrality";
+ filename += gCentrality; filename += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
+ else filename += "All.PttFrom";
+ }
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ 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
+ TFile *f = 0x0;
+ if(!gHistBalanceFunction2D) {
+ TFile::Open(filename.Data());
+ if((!f)||(!f->IsOpen())) {
+ Printf("The file %s is not found. Aborting...",filename);
+ return listBF;
+ }
+ //f->ls();
+ }
+
+ //Latex
+ TString centralityLatex = "Centrality: ";
+ centralityLatex += centralityArray[gCentrality-1];
+ centralityLatex += "%";
+
+ TString psiLatex;
+ if((psiMin == -0.5)&&(psiMax == 0.5))
+ psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
+ else if((psiMin == 0.5)&&(psiMax == 1.5))
+ psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
+ else if((psiMin == 1.5)&&(psiMax == 2.5))
+ psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
+ else
+ psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
+
+ TString pttLatex = Form("%.1f",ptTriggerMin);
+ pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
+ pttLatex += " GeV/c";
+
+ TString ptaLatex = Form("%.1f",ptAssociatedMin);
+ ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
+ ptaLatex += " GeV/c";
+
+ TLatex *latexInfo1 = new TLatex();
+ latexInfo1->SetNDC();
+ latexInfo1->SetTextSize(0.045);
+ latexInfo1->SetTextColor(1);
+
+
+ //============================================================//
+ //Get subtracted and mixed balance function
+ TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
+ TH2D *gHistBalanceFunctionMixed2D = 0x0;
+ if(!gHistBalanceFunction2D) {
+ gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
+ gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
+ }
+ else {
+ gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
+ gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
+ }
+
+ TH1D *gHistBalanceFunctionSubtracted = NULL;
+ TH1D *gHistBalanceFunctionMixed = NULL;
+
+ TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
+ TH1D *gHistBalanceFunctionMixed_scale = NULL;
+
+ if(kProjectInEta){
+ gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
+ gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
+ gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
+ 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("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
+ gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
+ gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
+ gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
+ gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
+ gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
+ }
+
+ gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
+ gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
+
+ gHistBalanceFunctionMixed->SetMarkerStyle(25);
+
+ TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
+ c1->SetFillColor(10);
+ c1->SetHighLightColor(10);
+ c1->SetLeftMargin(0.15);
+ gHistBalanceFunctionSubtracted->DrawCopy("E");
+ gHistBalanceFunctionMixed->DrawCopy("ESAME");
- latexInfo1->DrawLatex(0.68,0.88,centralityLatex.Data());
- latexInfo1->DrawLatex(0.68,0.82,psiLatex.Data());
- latexInfo1->DrawLatex(0.68,0.76,pttLatex.Data());
- latexInfo1->DrawLatex(0.68,0.70,ptaLatex.Data());
+ legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
+ legend->SetTextSize(0.045);
+ legend->SetTextFont(42);
+ legend->SetBorderSize(0);
+ legend->SetFillStyle(0);
+ legend->SetFillColor(10);
+ legend->SetMargin(0.25);
+ legend->SetShadowColor(10);
+ legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
+ legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
+ legend->Draw();
+
+
+ TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
+
+ if(bRootMoments){
+ meanLatex = "#mu = ";
+ meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
+ meanLatex += " #pm ";
+ meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
- TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
- c4->SetFillColor(10); c4->SetHighLightColor(10);
- c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
- gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
- gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
- gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
- gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
- gHistBalanceFunctionSubtracted->DrawCopy("lego2");
-
- latexInfo1->DrawLatex(0.68,0.88,centralityLatex.Data());
- latexInfo1->DrawLatex(0.68,0.82,psiLatex.Data());
- latexInfo1->DrawLatex(0.68,0.76,pttLatex.Data());
- latexInfo1->DrawLatex(0.68,0.70,ptaLatex.Data());
+ 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));
+
+
+ // store in txt files
+ TString meanFileName = filename;
+ if(kProjectInEta)
+ meanFileName= "deltaEtaProjection_Mean.txt";
+ //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
+ else
+ meanFileName = "deltaPhiProjection_Mean.txt";
+ //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
+ ofstream fileMean(meanFileName.Data(),ios::app);
+ fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
+ fileMean.close();
+
+ TString rmsFileName = filename;
+ if(kProjectInEta)
+ rmsFileName = "deltaEtaProjection_Rms.txt";
+ //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
+ else
+ rmsFileName = "deltaPhiProjection_Rms.txt");
+ //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
+ ofstream fileRms(rmsFileName.Data(),ios::app);
+ fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
+ fileRms.close();
+
+ TString skewnessFileName = filename;
+ if(kProjectInEta)
+ skewnessFileName = "deltaEtaProjection_Skewness.txt";
+ //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
+ else
+ skewnessFileName = "deltaPhiProjection_Skewness.txt";
+ //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
+ ofstream fileSkewness(skewnessFileName.Data(),ios::app);
+ fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
+ fileSkewness.close();
+
+ TString kurtosisFileName = filename;
+ if(kProjectInEta)
+ kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
+ //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
+ else
+ kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
+ //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
+ ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
+ fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
+ fileKurtosis.close();
+ }
+ // calculate the moments by hand
+ else{
+
+ Double_t meanAnalytical, meanAnalyticalError;
+ Double_t sigmaAnalytical, sigmaAnalyticalError;
+ Double_t skewnessAnalytical, skewnessAnalyticalError;
+ Double_t kurtosisAnalytical, kurtosisAnalyticalError;
+
+ Int_t gDeltaEtaPhi = 2;
+ if(kProjectInEta) gDeltaEtaPhi = 1;
+
+ AliBalancePsi *bHelper = new AliBalancePsi;
+ bHelper->GetMomentsAnalytical(gDeltaEtaPhi,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);
+
+ // store in txt files
+ TString meanFileName = filename;
+ if(kProjectInEta)
+ meanFileName = "deltaEtaProjection_Mean.txt";
+ //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
+ else
+ meanFileName = "deltaPhiProjection_Mean.txt";
+ //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
+ ofstream fileMean(meanFileName.Data(),ios::app);
+ fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
+ fileMean.close();
+
+ TString rmsFileName = filename;
+ if(kProjectInEta)
+ rmsFileName = "deltaEtaProjection_Rms.txt";
+ //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
+ else
+ rmsFileName = "deltaPhiProjection_Rms.txt";
+//rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
+ ofstream fileRms(rmsFileName.Data(),ios::app);
+ fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
+ fileRms.close();
+
+ TString skewnessFileName = filename;
+ if(kProjectInEta)
+ skewnessFileName = "deltaEtaProjection_Skewness.txt";
+ //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
+ else
+ skewnessFileName = "deltaPhiProjection_Skewness.txt";
+ //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
+ ofstream fileSkewness(skewnessFileName.Data(),ios::app);
+ fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
+ fileSkewness.close();
+
+ TString kurtosisFileName = filename;
+ if(kProjectInEta)
+ kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
+ //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
+ else
+ kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
+ //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
+ ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
+ fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
+ fileKurtosis.close();
+ }
+
+ // Weighted mean as calculated for 1D analysis
+ Double_t weightedMean, weightedMeanError;
+ GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError);
+ Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
+
+ // store in txt files
+ TString weightedMeanFileName = filename;
+ if(kProjectInEta)
+ weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
+ //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
+ else
+ weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
+ //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
+ ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
+ fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
+ fileWeightedMean.close();
+
+
+ TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
+ c2->SetFillColor(10);
+ c2->SetHighLightColor(10);
+ c2->SetLeftMargin(0.15);
+ gHistBalanceFunctionSubtracted->DrawCopy("E");
+
+ 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());
+
+ TString pngName = filename;
+ if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
+ else pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
+
+ c2->SaveAs(pngName.Data());
+
+ TString outFileName = filename;
+ if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
+ else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
+ TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
+ gHistBalanceFunctionSubtracted->Write();
+ gHistBalanceFunctionMixed->Write();
+ fProjection->Close();
+}
+
+//____________________________________________________________//
+void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
+ Int_t gTrainNumber = 64,
+ const char* gCentralityEstimator = "V0M",
+ Int_t gBit = 128,
+ const char* gEventPlaneEstimator = "VZERO",
+ Int_t gCentrality = 1,
+ Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+ Double_t vertexZMin = -10.,
+ Double_t vertexZMax = 10.,
+ Double_t ptTriggerMin = -1.,
+ Double_t ptTriggerMax = -1.,
+ Double_t ptAssociatedMin = -1.,
+ Double_t ptAssociatedMax = -1.,
+ TString eventClass = "Multiplicity"
+) {
+ //Macro that draws the BF distributions for each centrality bin
+ //for reaction plane dependent analysis
+ //Author: Panos.Christakoglou@nikhef.nl
+ TGaxis::SetMaxDigits(3);
+ gStyle->SetPalette(55,0);
+
+ //Get the input file
+ TString filename = lhcPeriod;
+ if(lhcPeriod != ""){
+ //filename += "/Train"; filename += gTrainNumber;
+ filename +="/PttFrom";
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ filename += Form("%.1f",ptAssociatedMin); filename += "To";
+ filename += Form("%.1f",ptAssociatedMax);
+ filename += "/correlationFunction.";
+ }
+ else{
+ filename += "correlationFunction.";
+ }
+ if(eventClass == "Centrality"){
+ filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
+ filename += ".PsiAll.PttFrom";
+ }
+ else if(eventClass == "Multiplicity"){
+ filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
+ filename += ".PsiAll.PttFrom";
+ }
+ else{ // "EventPlane" (default)
+ filename += "Centrality";
+ filename += gCentrality; filename += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
+ else filename += "All.PttFrom";
+ }
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ filename += Form("%.1f",ptAssociatedMin); filename += "To";
+ filename += Form("%.1f",ptAssociatedMax);
+ filename += "_";
+ filename += Form("%.1f",psiMin);
+ filename += "-";
+ filename += Form("%.1f",psiMax);
+ filename += ".root";
+
+ //Open the file
+ TFile *f = TFile::Open(filename.Data());
+ if((!f)||(!f->IsOpen())) {
+ Printf("The file %s is not found. Aborting...",filename);
+ return listBF;
+ }
+ //f->ls();
+
+ TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
+ if(!gHistPN) return;
+ TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
+ if(!gHistNP) return;
+ TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
+ if(!gHistPP) return;
+ TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
+ if(!gHistNN) return;
+
+ gHistPN->Sumw2();
+ gHistPP->Sumw2();
+ gHistPN->Add(gHistPP,-1);
+ gHistNP->Sumw2();
+ gHistNN->Sumw2();
+ gHistNP->Add(gHistNN,-1);
+ gHistPN->Add(gHistNP);
+ gHistPN->Scale(0.5);
+ TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
+ gHistBalanceFunction2D->SetStats(kFALSE);
+ gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
+
+ //Draw the results
+ TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
+ c0->SetFillColor(10); c0->SetHighLightColor(10);
+ c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
+ gHistBalanceFunction2D->SetTitle("");
+ gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
+ gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
+ gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
+ gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
+ gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
+ gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
+ gHistBalanceFunction2D->DrawCopy("lego2");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
+
+ TString pttLatex = Form("%.1f",ptTriggerMin);
+ pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
+ pttLatex += " GeV/c";
+
+ TString ptaLatex = Form("%.1f",ptAssociatedMin);
+ ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
+ ptaLatex += " GeV/c";
+
+ TLatex *latexInfo1 = new TLatex();
+ latexInfo1->SetNDC();
+ latexInfo1->SetTextSize(0.045);
+ latexInfo1->SetTextColor(1);
+ latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
+ latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
+ latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
+
+ TString pngName = "BalanceFunction2D.";
+ pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
+ pngName += ".PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ pngName += ".png";
+
+ c0->SaveAs(pngName.Data());
+
+ drawProjections(gHistBalanceFunction2D,
+ kTRUE,
+ 1,36,
+ gCentrality,
+ psiMin,psiMax,
+ ptTriggerMin,ptTriggerMax,
+ ptAssociatedMin,ptAssociatedMax,
+ kTRUE,
+ eventClass,
+ kFALSE);
+
+ drawProjections(gHistBalanceFunction2D,
+ kFALSE,
+ 1,80,
+ gCentrality,
+ psiMin,psiMax,
+ ptTriggerMin,ptTriggerMax,
+ ptAssociatedMin,ptAssociatedMax,
+ kTRUE,
+ eventClass.Data(),
+ kFALSE);
+
+ TString outFileName = filename;
+ outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
+ gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
+ TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
+ gHistBalanceFunction2D->Write();
+ fOut->Close();
+
}
+
+//____________________________________________________________________//
+void GetWeightedMean1D(TH1D *gHistBalance, Bool_t kProjectInEta = kTRUE, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME) {
+
+ //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();
+ if(fStopBin > -1) fNumberOfBins = fStopBin;
+ Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
+ Double_t currentBinCenter = 0.;
+
+ for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
+
+ // in order to recover the |abs| in the 1D analysis
+ currentBinCenter = gHistBalance->GetBinCenter(i);
+ if(kProjectInEta){
+ if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
+ }
+ else{
+ if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
+ if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
+ }
+
+ gSumXi += currentBinCenter;
+ gSumBi += gHistBalance->GetBinContent(i);
+ gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
+ gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
+ gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
+ gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
+ gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,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)*(currentBinCenter*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) );
+
+ WM = delta;
+ WME = deltaError;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff