PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

   1 const Int_t numberOfCentralityBins = 9;
   2 TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
   3 Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.};
   4 Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.};
   5 TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
   6
   7 const Int_t gRebin = 1;
   8 void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
   9                               Double_t psiMin = 0., Double_t psiMax = 7.5,
  10                               Double_t ptTriggerMin = -1.,
  11                               Double_t ptTriggerMax = -1.,
  12                               Double_t ptAssociatedMin = -1.,
  13                               Double_t ptAssociatedMax = -1.) {
  14   //Macro that draws the BF distributions for each centrality bin
  15   //for reaction plane dependent analysis
  16   //Author: Panos.Christakoglou@nikhef.nl
  17   //Load the PWG2ebye library
  18   gSystem->Load("libANALYSIS.so");
  19   gSystem->Load("libANALYSISalice.so");
  20   gSystem->Load("libEventMixing.so");
  21   gSystem->Load("libCORRFW.so");
  22   gSystem->Load("libPWGTools.so");
  23   gSystem->Load("libPWGCFebye.so");
  24
  25   //Prepare the objects and return them
  26   TList *listBF = GetListOfObjects(filename);
  27   TList *listBFShuffled = GetListOfObjects(filename,kTRUE);
  28   if(!listBF) {
  29     Printf("The TList object was not created");
  30     return;
  31   }
  32   else
  33     draw(listBF,listBFShuffled,psiMin,psiMax,
  34          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
  35 }
  36
  37 //______________________________________________________//
  38 TList *GetListOfObjects(const char* filename, Bool_t kShuffling = kFALSE) {
  39   //Get the TList objects (QA, bf, bf shuffled)
  40   TList *listQA = 0x0;
  41   TList *listBF = 0x0;
  42   TList *listBFShuffling = 0x0;
  43
  44   //Open the file
  45   TFile *f = TFile::Open(filename);
  46   if((!f)||(!f->IsOpen())) {
  47     Printf("The file %s is not found. Aborting...",filename);
  48     return listBF;
  49   }
  50   //f->ls();
  51
  52   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  53   if(!dir) {
  54     Printf("The TDirectoryFile is not found. Aborting...",filename);
  55     return listBF;
  56   }
  57   //dir->ls();
  58
  59   TString listBFName;
  60   if(!kShuffling)
  61     listBFName = "listBFPsi";
  62   else if(kShuffling)
  63     listBFName = "listBFPsiShuffled";
  64   listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
  65   //listBF->ls();
  66
  67   //Get the histograms
  68   TString histoName;
  69   if(!kShuffling)
  70     histoName = "fHistPV0M";
  71   else if(kShuffling)
  72     histoName = "fHistP_shuffleV0M";
  73   AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  74   if(!fHistP) {
  75     Printf("fHistP %s not found!!!",histoName.Data());
  76     break;
  77   }
  78   fHistP->FillParent(); fHistP->DeleteContainers();
  79
  80   if(!kShuffling)
  81     histoName = "fHistNV0M";
  82   else if(kShuffling)
  83     histoName = "fHistN_shuffleV0M";
  84   AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  85   if(!fHistN) {
  86     Printf("fHistN %s not found!!!",histoName.Data());
  87     break;
  88   }
  89   fHistN->FillParent(); fHistN->DeleteContainers();
  90
  91   if(!kShuffling)
  92     histoName = "fHistPNV0M";
  93   else if(kShuffling)
  94     histoName = "fHistPN_shuffleV0M";
  95   AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  96   if(!fHistPN) {
  97     Printf("fHistPN %s not found!!!",histoName.Data());
  98     break;
  99   }
 100   fHistPN->FillParent(); fHistPN->DeleteContainers();
 101
 102   if(!kShuffling)
 103     histoName = "fHistNPV0M";
 104   else if(kShuffling)
 105     histoName = "fHistNP_shuffleV0M";
 106   AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 107   if(!fHistNP) {
 108     Printf("fHistNP %s not found!!!",histoName.Data());
 109     break;
 110   }
 111   fHistNP->FillParent(); fHistNP->DeleteContainers();
 112
 113   if(!kShuffling)
 114     histoName = "fHistPPV0M";
 115   else if(kShuffling)
 116     histoName = "fHistPP_shuffleV0M";
 117   AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 118   if(!fHistPP) {
 119     Printf("fHistPP %s not found!!!",histoName.Data());
 120     break;
 121   }
 122   fHistPP->FillParent(); fHistPP->DeleteContainers();
 123
 124   if(!kShuffling)
 125     histoName = "fHistNNV0M";
 126   else if(kShuffling)
 127     histoName = "fHistNN_shuffleV0M";
 128   AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 129   if(!fHistNN) {
 130     Printf("fHistNN %s not found!!!",histoName.Data());
 131     break;
 132   }
 133   fHistNN->FillParent(); fHistNN->DeleteContainers();
 134
 135   return listBF;
 136 }
 137
 138 //______________________________________________________//
 139 void draw(TList *listBF, TList *listBFShuffled,
 140           Double_t psiMin, Double_t psiMax,
 141           Double_t ptTriggerMin, Double_t ptTriggerMax,
 142           Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
 143   gROOT->LoadMacro("~/SetPlotStyle.C");
 144   SetPlotStyle();
 145   gStyle->SetPalette(1,0);
 146
 147   //balance function
 148   AliTHn *hP = NULL;
 149   AliTHn *hN = NULL;
 150   AliTHn *hPN = NULL;
 151   AliTHn *hNP = NULL;
 152   AliTHn *hPP = NULL;
 153   AliTHn *hNN = NULL;
 154   //listBF->ls();
 155   //Printf("=================");
 156   hP = (AliTHn*) listBF->FindObject("fHistPV0M");
 157   hN = (AliTHn*) listBF->FindObject("fHistNV0M");
 158   hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
 159   hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
 160   hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
 161   hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
 162
 163   AliBalancePsi *b = new AliBalancePsi();
 164   b->SetHistNp(hP);
 165   b->SetHistNn(hN);
 166   b->SetHistNpn(hPN);
 167   b->SetHistNnp(hNP);
 168   b->SetHistNpp(hPP);
 169   b->SetHistNnn(hNN);
 170
 171   //balance function shuffling
 172   AliTHn *hPShuffled = NULL;
 173   AliTHn *hNShuffled = NULL;
 174   AliTHn *hPNShuffled = NULL;
 175   AliTHn *hNPShuffled = NULL;
 176   AliTHn *hPPShuffled = NULL;
 177   AliTHn *hNNShuffled = NULL;
 178   //listBFShuffled->ls();
 179   hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
 180   hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
 181   hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
 182   hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
 183   hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
 184   hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
 185
 186   AliBalancePsi *bShuffled = new AliBalancePsi();
 187   bShuffled->SetHistNp(hPShuffled);
 188   bShuffled->SetHistNn(hNShuffled);
 189   bShuffled->SetHistNpn(hPNShuffled);
 190   bShuffled->SetHistNnp(hNPShuffled);
 191   bShuffled->SetHistNpp(hPPShuffled);
 192   bShuffled->SetHistNnn(hNNShuffled);
 193
 194   TH2D *gHistBalanceFunction;
 195   TH2D *gHistBalanceFunctionShuffled;
 196   TCanvas *c1;
 197   TString histoTitle, pngName;
 198   TLegend *legend;
 199
 200   //loop over the centrality bins
 201   //for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) {
 202   histoTitle = "Centrality: ";
 203   histoTitle += centralityArray[6];
 204   histoTitle += "%";
 205   histoTitle += " | "; histoTitle += psiMin;
 206   histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
 207
 208   gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 209   gHistBalanceFunction->SetTitle(histoTitle.Data());
 210   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 211
 212   gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 213
 214   c1 = new TCanvas(histoTitle.Data(),"",0,0,600,500);
 215   c1->SetFillColor(10);
 216   c1->SetHighLightColor(10);
 217   c1->SetLeftMargin(0.15);
 218   gHistBalanceFunction->Draw("lego2");
 219   //gHistBalanceFunctionShuffled->Draw("ESAME");
 220
 221   legend = new TLegend(0.18,0.6,0.45,0.82,"","brNDC");
 222   legend->SetTextSize(0.045);
 223   legend->SetTextFont(42);
 224   legend->SetBorderSize(0);
 225   legend->SetFillStyle(0);
 226   legend->SetFillColor(10);
 227   legend->SetMargin(0.25);
 228   legend->SetShadowColor(10);
 229   legend->AddEntry(gHistBalanceFunction,"Data","lp");
 230   legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
 231   legend->Draw();
 232
 233   pngName = "BalanceFunctionDeltaEta.Centrality";
 234   pngName += centralityArray[6];
 235   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 236   pngName += ".png";
 237   c1->SaveAs(pngName.Data());
 238 }
 239
 240 //____________________________________________________________________//
 241 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
 242   //Prints the calculated width of the BF and its error
 243   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 244   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 245   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 246   Double_t deltaBalP2 = 0.0, integral = 0.0;
 247   Double_t deltaErrorNew = 0.0;
 248
 249   //Retrieve this variables from Histogram
 250   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
 251   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
 252
 253   cout<<"=================================================="<<endl;
 254   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
 255   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
 256   cout<<"=================================================="<<endl;
 257   for(Int_t i = 1; i <= fNumberOfBins; i++) {
 258
 259     // this is to simulate |Delta eta| or |Delta phi|
 260     if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
 261
 262     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
 263
 264     gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
 265     gSumBi += gHistBalance->GetBinContent(i);
 266     gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
 267     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 268     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 269     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 270     gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 271
 272     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
 273     integral += fP2Step*gHistBalance->GetBinContent(i);
 274   }
 275   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
 276     deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 277
 278   Double_t integralError = TMath::Sqrt(deltaBalP2);
 279
 280   Double_t delta = gSumBiXi / gSumBi;
 281   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 282   cout<<"=================================================="<<endl;
 283   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 284   cout<<"New error: "<<deltaErrorNew<<endl;
 285   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 286   cout<<"=================================================="<<endl;
 287   cout<<endl;
 288 }