PWGCF/EBYE/macros/drawBalanceFunctionPsi.C

   1 const Int_t numberOfCentralityBins = 10;
   2 TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-1","1-2"};
   3
   4 void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root",
   5                             Int_t gCentrality = 1,
   6                             Int_t gDeltaEtaDeltaPhi = 2,
   7                             Double_t psiMin = -0.5, Double_t psiMax = 0.5,
   8                             Double_t ptTriggerMin = -1.,
   9                             Double_t ptTriggerMax = -1.,
  10                             Double_t ptAssociatedMin = -1.,
  11                             Double_t ptAssociatedMax = -1.) {
  12   //Macro that draws the BF distributions for each centrality bin
  13   //for reaction plane dependent analysis
  14   //Author: Panos.Christakoglou@nikhef.nl
  15   //Load the PWG2ebye library
  16   gSystem->Load("libANALYSIS.so");
  17   gSystem->Load("libANALYSISalice.so");
  18   gSystem->Load("libEventMixing.so");
  19   gSystem->Load("libCORRFW.so");
  20   gSystem->Load("libPWGTools.so");
  21   gSystem->Load("libPWGCFebye.so");
  22
  23   //Prepare the objects and return them
  24   TList *listBF = GetListOfObjects(filename,gCentrality,0);
  25   TList *listBFShuffled = GetListOfObjects(filename,gCentrality,1);
  26   TList *listBFMixed = GetListOfObjects(filename,gCentrality,2);
  27   if(!listBF) {
  28     Printf("The TList object was not created");
  29     return;
  30   }
  31   else
  32     draw(listBF,listBFShuffled,listBFMixed,
  33          gCentrality,gDeltaEtaDeltaPhi,
  34          psiMin,psiMax,
  35          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
  36 }
  37
  38 //______________________________________________________//
  39 TList *GetListOfObjects(const char* filename,
  40                         Int_t gCentrality, Int_t kData = 1) {
  41   //Get the TList objects (QA, bf, bf shuffled)
  42   //kData == 0: data
  43   //kData == 1: shuffling
  44   //kData == 2: mixing
  45   TList *listQA = 0x0;
  46   TList *listBF = 0x0;
  47
  48   //Open the file
  49   TFile *f = TFile::Open(filename);
  50   if((!f)||(!f->IsOpen())) {
  51     Printf("The file %s is not found. Aborting...",filename);
  52     return listBF;
  53   }
  54   f->ls();
  55
  56   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  57   if(!dir) {
  58     Printf("The TDirectoryFile is not found. Aborting...",filename);
  59     return listBF;
  60   }
  61   dir->ls();
  62
  63   TString listBFName;
  64   if(kData == 0) {
  65     //cout<<"no shuffling - no mixing"<<endl;
  66     listBFName = "listBFPsi_";
  67   }
  68   else if(kData == 1) {
  69     //cout<<"shuffling - no mixing"<<endl;
  70     listBFName = "listBFPsiShuffled_";
  71   }
  72   else if(kData == 2) {
  73     //cout<<"no shuffling - mixing"<<endl;
  74     listBFName = "listBFPsiMixed_";
  75   }
  76   listBFName += centralityArray[gCentrality-1];
  77   listBFName += "_Bit128_V0M";
  78   listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
  79   cout<<"======================================================="<<endl;
  80   cout<<"List name: "<<listBF->GetName()<<endl;
  81   listBF->ls();
  82
  83   //Get the histograms
  84   TString histoName;
  85   if(kData == 0)
  86     histoName = "fHistPV0M";
  87   else if(kData == 1)
  88     histoName = "fHistP_shuffleV0M";
  89   else if(kData == 2)
  90     histoName = "fHistPV0M";
  91   AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  92   if(!fHistP) {
  93     Printf("fHistP %s not found!!!",histoName.Data());
  94     break;
  95   }
  96   fHistP->FillParent(); fHistP->DeleteContainers();
  97
  98   if(kData == 0)
  99     histoName = "fHistNV0M";
 100   if(kData == 1)
 101     histoName = "fHistN_shuffleV0M";
 102   if(kData == 2)
 103     histoName = "fHistNV0M";
 104   AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 105   if(!fHistN) {
 106     Printf("fHistN %s not found!!!",histoName.Data());
 107     break;
 108   }
 109   fHistN->FillParent(); fHistN->DeleteContainers();
 110
 111   if(kData == 0)
 112     histoName = "fHistPNV0M";
 113   if(kData == 1)
 114     histoName = "fHistPN_shuffleV0M";
 115   if(kData == 2)
 116     histoName = "fHistPNV0M";
 117   AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 118   if(!fHistPN) {
 119     Printf("fHistPN %s not found!!!",histoName.Data());
 120     break;
 121   }
 122   fHistPN->FillParent(); fHistPN->DeleteContainers();
 123
 124   if(kData == 0)
 125     histoName = "fHistNPV0M";
 126   if(kData == 1)
 127     histoName = "fHistNP_shuffleV0M";
 128   if(kData == 2)
 129     histoName = "fHistNPV0M";
 130   AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 131   if(!fHistNP) {
 132     Printf("fHistNP %s not found!!!",histoName.Data());
 133     break;
 134   }
 135   fHistNP->FillParent(); fHistNP->DeleteContainers();
 136
 137   if(kData == 0)
 138     histoName = "fHistPPV0M";
 139   if(kData == 1)
 140     histoName = "fHistPP_shuffleV0M";
 141   if(kData == 2)
 142     histoName = "fHistPPV0M";
 143   AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 144   if(!fHistPP) {
 145     Printf("fHistPP %s not found!!!",histoName.Data());
 146     break;
 147   }
 148   fHistPP->FillParent(); fHistPP->DeleteContainers();
 149
 150   if(kData == 0)
 151     histoName = "fHistNNV0M";
 152   if(kData == 1)
 153     histoName = "fHistNN_shuffleV0M";
 154   if(kData == 2)
 155     histoName = "fHistNNV0M";
 156   AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 157   if(!fHistNN) {
 158     Printf("fHistNN %s not found!!!",histoName.Data());
 159     break;
 160   }
 161   fHistNN->FillParent(); fHistNN->DeleteContainers();
 162
 163   return listBF;
 164 }
 165
 166 //______________________________________________________//
 167 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 168           Int_t gCentrality, Int_t gDeltaEtaDeltaPhi,
 169           Double_t psiMin, Double_t psiMax,
 170           Double_t ptTriggerMin, Double_t ptTriggerMax,
 171           Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
 172   gROOT->LoadMacro("~/SetPlotStyle.C");
 173   SetPlotStyle();
 174   gStyle->SetPalette(1,0);
 175
 176   const Int_t gRebin = gDeltaEtaDeltaPhi; //rebin by 2 the Delta phi projection
 177
 178   //balance function
 179   AliTHn *hP = NULL;
 180   AliTHn *hN = NULL;
 181   AliTHn *hPN = NULL;
 182   AliTHn *hNP = NULL;
 183   AliTHn *hPP = NULL;
 184   AliTHn *hNN = NULL;
 185   //listBF->ls();
 186   //Printf("=================");
 187   hP = (AliTHn*) listBF->FindObject("fHistPV0M");
 188   hN = (AliTHn*) listBF->FindObject("fHistNV0M");
 189   hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
 190   hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
 191   hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
 192   hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
 193
 194   AliBalancePsi *b = new AliBalancePsi();
 195   b->SetHistNp(hP);
 196   b->SetHistNn(hN);
 197   b->SetHistNpn(hPN);
 198   b->SetHistNnp(hNP);
 199   b->SetHistNpp(hPP);
 200   b->SetHistNnn(hNN);
 201
 202   //balance function shuffling
 203   AliTHn *hPShuffled = NULL;
 204   AliTHn *hNShuffled = NULL;
 205   AliTHn *hPNShuffled = NULL;
 206   AliTHn *hNPShuffled = NULL;
 207   AliTHn *hPPShuffled = NULL;
 208   AliTHn *hNNShuffled = NULL;
 209   //listBFShuffled->ls();
 210   hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
 211   hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
 212   hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
 213   hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
 214   hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
 215   hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
 216
 217   AliBalancePsi *bShuffled = new AliBalancePsi();
 218   bShuffled->SetHistNp(hPShuffled);
 219   bShuffled->SetHistNn(hNShuffled);
 220   bShuffled->SetHistNpn(hPNShuffled);
 221   bShuffled->SetHistNnp(hNPShuffled);
 222   bShuffled->SetHistNpp(hPPShuffled);
 223   bShuffled->SetHistNnn(hNNShuffled);
 224
 225   //balance function mixing
 226   AliTHn *hPMixed = NULL;
 227   AliTHn *hNMixed = NULL;
 228   AliTHn *hPNMixed = NULL;
 229   AliTHn *hNPMixed = NULL;
 230   AliTHn *hPPMixed = NULL;
 231   AliTHn *hNNMixed = NULL;
 232   //listBFMixed->ls();
 233   hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
 234   hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
 235   hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
 236   hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
 237   hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
 238   hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
 239
 240   AliBalancePsi *bMixed = new AliBalancePsi();
 241   bMixed->SetHistNp(hPMixed);
 242   bMixed->SetHistNn(hNMixed);
 243   bMixed->SetHistNpn(hPNMixed);
 244   bMixed->SetHistNnp(hNPMixed);
 245   bMixed->SetHistNpp(hPPMixed);
 246   bMixed->SetHistNnn(hNNMixed);
 247
 248   TH1D *gHistBalanceFunction;
 249   TH1D *gHistBalanceFunctionShuffled;
 250   TH1D *gHistBalanceFunctionMixed;
 251   TH1D *gHistBalanceFunctionSubtracted;
 252   TString histoTitle, pngName;
 253   TLegend *legend;
 254
 255   histoTitle = "Centrality: ";
 256   histoTitle += centralityArray[gCentrality-1];
 257   histoTitle += "%";
 258   if((psiMin == -0.5)&&(psiMax == 0.5))
 259     histoTitle += " (-7.5^{o} < #phi - #Psi_{2} < 7.5^{o})";
 260   else if((psiMin == 0.5)&&(psiMax == 1.5))
 261     histoTitle += " (37.5^{o} < #phi - #Psi_{2} < 52.5^{o})";
 262   else if((psiMin == 1.5)&&(psiMax == 2.5))
 263     histoTitle += " (82.5^{o} < #phi - #Psi_{2} < 97.5^{o})";
 264   else
 265     histoTitle += " (0^{o} < #phi - #Psi_{2} < 180^{o})";
 266
 267   //Raw balance function
 268   gHistBalanceFunction = b->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 269   gHistBalanceFunction->SetMarkerStyle(20);
 270   gHistBalanceFunction->SetTitle(histoTitle.Data());
 271   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 272   gHistBalanceFunction->SetName("gHistBalanceFunction");
 273
 274   //Shuffled balance function
 275   gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 276   gHistBalanceFunctionShuffled->SetMarkerStyle(24);
 277   gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 278
 279   //Mixed balance function
 280   gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 281   gHistBalanceFunctionMixed->SetMarkerStyle(25);
 282   gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 283
 284   //Subtracted balance function
 285   gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunction->Clone());
 286   gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 287   gHistBalanceFunctionSubtracted->Rebin(gRebin);
 288   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
 289   gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 290   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 291   gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 292
 293   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 294   c1->SetFillColor(10);
 295   c1->SetHighLightColor(10);
 296   c1->SetLeftMargin(0.15);
 297   gHistBalanceFunction->DrawCopy("E");
 298   gHistBalanceFunctionShuffled->DrawCopy("ESAME");
 299   gHistBalanceFunctionMixed->DrawCopy("ESAME");
 300
 301   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
 302   legend->SetTextSize(0.045);
 303   legend->SetTextFont(42);
 304   legend->SetBorderSize(0);
 305   legend->SetFillStyle(0);
 306   legend->SetFillColor(10);
 307   legend->SetMargin(0.25);
 308   legend->SetShadowColor(10);
 309   legend->AddEntry(gHistBalanceFunction,"Data","lp");
 310   legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
 311   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
 312   legend->Draw();
 313
 314   pngName = "BalanceFunctionDeltaEta.Centrality";
 315   pngName += centralityArray[gCentrality-1];
 316   pngName += ".Psi"; //pngName += psiMin; pngName += "To"; pngName += psiMax;
 317   pngName += ".png";
 318   c1->SaveAs(pngName.Data());
 319
 320   GetWeightedMean(gHistBalanceFunction);
 321   GetWeightedMean(gHistBalanceFunctionShuffled);
 322
 323   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
 324   meanLatex = "#mu = ";
 325   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
 326   meanLatex += " #pm ";
 327   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
 328
 329   rmsLatex = "#sigma = ";
 330   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
 331   rmsLatex += " #pm ";
 332   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
 333
 334   skewnessLatex = "S = ";
 335   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
 336   skewnessLatex += " #pm ";
 337   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
 338
 339   kurtosisLatex = "K = ";
 340   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
 341   kurtosisLatex += " #pm ";
 342   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
 343   Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
 344   Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
 345   Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
 346   Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
 347
 348   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
 349   c2->SetFillColor(10);
 350   c2->SetHighLightColor(10);
 351   c2->SetLeftMargin(0.15);
 352   gHistBalanceFunctionSubtracted->DrawCopy("E");
 353
 354   TLatex *latex = new TLatex();
 355   latex->SetNDC();
 356   latex->SetTextSize(0.035);
 357   latex->SetTextColor(1);
 358   latex->DrawLatex(0.64,0.85,meanLatex.Data());
 359   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
 360   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
 361   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
 362
 363   TString newFileName = "balanceFunction.Centrality";
 364   newFileName += gCentrality; newFileName += ".In";
 365   if(gDeltaEtaDeltaPhi == 1) newFileName += "DeltaEta.Psi";
 366   else if(gDeltaEtaDeltaPhi == 2) newFileName += "DeltaPhi.Psi";
 367   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 368   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 369   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 370   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 371   else newFileName += "All.PttFrom";
 372   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 373   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 374   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 375   newFileName += Form("%.1f",ptAssociatedMax);
 376   newFileName += ".root";
 377
 378   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 379   fOutput->cd();
 380   gHistBalanceFunction->Write();
 381   gHistBalanceFunctionShuffled->Write();
 382   gHistBalanceFunctionMixed->Write();
 383   gHistBalanceFunctionSubtracted->Write();
 384   fOutput->Close();
 385 }
 386
 387 //____________________________________________________________________//
 388 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
 389   //Prints the calculated width of the BF and its error
 390   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 391   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 392   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 393   Double_t deltaBalP2 = 0.0, integral = 0.0;
 394   Double_t deltaErrorNew = 0.0;
 395
 396   //Retrieve this variables from Histogram
 397   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
 398   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
 399
 400   cout<<"=================================================="<<endl;
 401   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
 402   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
 403   cout<<"=================================================="<<endl;
 404   for(Int_t i = 1; i <= fNumberOfBins; i++) {
 405     // this is to simulate |Delta eta| or |Delta phi|
 406     if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
 407
 408     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
 409
 410     gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
 411     gSumBi += gHistBalance->GetBinContent(i);
 412     gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
 413     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 414     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 415     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 416     gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 417
 418     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
 419     integral += fP2Step*gHistBalance->GetBinContent(i);
 420   }
 421   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
 422     deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 423
 424   Double_t integralError = TMath::Sqrt(deltaBalP2);
 425
 426   Double_t delta = gSumBiXi / gSumBi;
 427   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 428   cout<<"=================================================="<<endl;
 429   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 430   cout<<"New error: "<<deltaErrorNew<<endl;
 431   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 432   cout<<"=================================================="<<endl;
 433   cout<<endl;
 434 }
 435
 436 //______________________________________________________//
 437 void drawBFPsi(const char* lhcPeriod = "LHC10h",
 438                Int_t gCentrality = 1,
 439                Int_t gDeltaEtaDeltaPhi = 2,
 440                Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 441                Double_t ptTriggerMin = -1.,
 442                Double_t ptTriggerMax = -1.,
 443                Double_t ptAssociatedMin = -1.,
 444                Double_t ptAssociatedMax = -1.) {
 445   //Macro that draws the BF distributions for each centrality bin
 446   //for reaction plane dependent analysis
 447   //Author: Panos.Christakoglou@nikhef.nl
 448   gROOT->LoadMacro("~/SetPlotStyle.C");
 449   SetPlotStyle();
 450
 451   //Get the input file
 452   TString filename = lhcPeriod; filename +="/PttFrom";
 453   filename += Form("%.1f",ptTriggerMin); filename += "To";
 454   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 455   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 456   filename += Form("%.1f",ptAssociatedMax);
 457   filename += "/balanceFunction.Centrality";
 458   filename += gCentrality; filename += ".In";
 459   if(gDeltaEtaDeltaPhi == 1) filename += "DeltaEta.Psi";
 460   else if(gDeltaEtaDeltaPhi == 2) filename += "DeltaPhi.Psi";
 461   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 462   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 463   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 464   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 465   else filename += "All.PttFrom";
 466   filename += Form("%.1f",ptTriggerMin); filename += "To";
 467   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 468   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 469   filename += Form("%.1f",ptAssociatedMax);  filename += ".root";
 470
 471   //Open the file
 472   TFile *f = TFile::Open(filename.Data());
 473   if((!f)||(!f->IsOpen())) {
 474     Printf("The file %s is not found. Aborting...",filename);
 475     return listBF;
 476   }
 477   //f->ls();
 478
 479   //Raw balance function
 480   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 481   gHistBalanceFunction->SetStats(kFALSE);
 482   gHistBalanceFunction->SetMarkerStyle(20);
 483   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 484   if(gDeltaEtaDeltaPhi == 2) {
 485     gHistBalanceFunction->GetYaxis()->SetTitle("B(#Delta #varphi)");
 486     gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
 487   }
 488
 489   //Shuffled balance function
 490   TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 491   gHistBalanceFunction->SetStats(kFALSE);
 492   gHistBalanceFunctionShuffled->SetMarkerStyle(24);
 493
 494   //Mixed balance function
 495   TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 496   gHistBalanceFunction->SetStats(kFALSE);
 497   gHistBalanceFunctionMixed->SetMarkerStyle(25);
 498
 499   //Subtracted balance function
 500   TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 501   gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 502   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
 503   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 504   if(gDeltaEtaDeltaPhi == 2) {
 505     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("B(#Delta #varphi)");
 506     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
 507   }
 508
 509   TString pngName;
 510   TLegend *legend;
 511
 512   TString centralityLatex = "Centrality: ";
 513   centralityLatex += centralityArray[gCentrality-1];
 514   centralityLatex += "%";
 515
 516   TString psiLatex;
 517   if((psiMin == -0.5)&&(psiMax == 0.5))
 518     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 519   else if((psiMin == 0.5)&&(psiMax == 1.5))
 520     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 521   else if((psiMin == 1.5)&&(psiMax == 2.5))
 522     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 523   else
 524     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 525
 526   TString pttLatex = Form("%.1f",ptTriggerMin);
 527   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 528   pttLatex += " GeV/c";
 529
 530   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 531   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 532   ptaLatex += " GeV/c";
 533
 534   //Draw the results
 535   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 536   c1->SetFillColor(10); c1->SetHighLightColor(10);
 537   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 538   gHistBalanceFunction->SetTitle("");
 539   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 540   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 541   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 542   gHistBalanceFunction->DrawCopy("E");
 543   gHistBalanceFunctionShuffled->DrawCopy("ESAME");
 544   gHistBalanceFunctionMixed->DrawCopy("ESAME");
 545
 546   legend = new TLegend(0.2,0.72,0.45,0.92,"","brNDC");
 547   legend->SetTextSize(0.05);
 548   legend->SetTextFont(42);
 549   legend->SetBorderSize(0);
 550   legend->SetFillStyle(0);
 551   legend->SetFillColor(10);
 552   legend->SetMargin(0.25);
 553   legend->SetShadowColor(10);
 554   legend->AddEntry(gHistBalanceFunction,"Data","lp");
 555   legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
 556   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
 557   legend->Draw();
 558
 559   TLatex *latexInfo1 = new TLatex();
 560   latexInfo1->SetNDC();
 561   latexInfo1->SetTextSize(0.04);
 562   latexInfo1->SetTextColor(1);
 563   latexInfo1->DrawLatex(0.58,0.88,centralityLatex.Data());
 564   latexInfo1->DrawLatex(0.58,0.82,psiLatex.Data());
 565   latexInfo1->DrawLatex(0.58,0.76,pttLatex.Data());
 566   latexInfo1->DrawLatex(0.58,0.70,ptaLatex.Data());
 567
 568   //pngName = "BalanceFunctionDeltaEta.Centrality";
 569   //pngName += centralityArray[gCentrality-1];
 570   //pngName += ".Psi"; //pngName += psiMin; pngName += "To"; pngName += psiMax;
 571   //pngName += ".png";
 572   //c1->SaveAs(pngName.Data());
 573
 574   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
 575   c2->SetFillColor(10); c2->SetHighLightColor(10);
 576   c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 577   gHistBalanceFunctionSubtracted->SetTitle("");
 578   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 579   gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 580   gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 581   gHistBalanceFunctionSubtracted->DrawCopy("E");
 582
 583   //Opening the output ascii files
 584   TString filenameMean = "meanIn";
 585   TString filenameSigma = "sigmaIn";
 586   TString filenameSkewness = "skewnessIn";
 587   TString filenameKurtosis = "kurtosisIn";
 588   if(gDeltaEtaDeltaPhi == 1) {
 589     filenameMean += "DeltaEta.Psi"; filenameSigma += "DeltaEta.Psi";
 590     filenameSkewness += "DeltaEta.Psi"; filenameKurtosis += "DeltaEta.Psi";
 591   }
 592   else if(gDeltaEtaDeltaPhi == 2) {
 593     filenameMean += "DeltaPhi.Psi"; filenameSigma += "DeltaPhi.Psi";
 594     filenameSkewness += "DeltaPhi.Psi"; filenameKurtosis += "DeltaPhi.Psi";
 595   }
 596   if((psiMin == -0.5)&&(psiMax == 0.5)) {
 597     filenameMean += "InPlane.Ptt"; filenameSigma += "InPlane.Ptt";
 598     filenameSkewness += "InPlane.Ptt"; filenameKurtosis += "InPlane.Ptt";
 599   }
 600   else if((psiMin == 0.5)&&(psiMax == 1.5)) {
 601     filenameMean += "Intermediate.Ptt"; filenameSigma += "Intermediate.Ptt";
 602     filenameSkewness += "Intermediate.Ptt";
 603     filenameKurtosis += "Intermediate.Ptt";
 604   }
 605   else if((psiMin == 1.5)&&(psiMax == 2.5)) {
 606     filenameMean += "OutOfPlane.Ptt"; filenameSigma += "OutOfPlane.Ptt";
 607     filenameSkewness += "OutOfPlane.Ptt";
 608     filenameKurtosis += "OutOfPlane.Ptt";
 609   }
 610   else if((psiMin == 2.5)&&(psiMax == 3.5)) {
 611     filenameMean += "Rest.Ptt"; filenameSigma += "Rest.Ptt";
 612     filenameSkewness += "Rest.Ptt"; filenameKurtosis += "Rest.Ptt";
 613   }
 614   else {
 615     filenameMean += "All.Ptt"; filenameSigma += "All.Ptt";
 616     filenameSkewness += "All.Ptt"; filenameKurtosis += "All.Ptt";
 617   }
 618   filenameMean += Form("%.1f",ptTriggerMin); filenameMean += "To";
 619   filenameMean += Form("%.1f",ptTriggerMax); filenameMean += "PtaFrom";
 620   filenameMean += Form("%.1f",ptAssociatedMin); filenameMean += "To";
 621   filenameMean += Form("%.1f",ptAssociatedMax);  filenameMean += ".txt";
 622   filenameSigma += Form("%.1f",ptTriggerMin); filenameSigma += "To";
 623   filenameSigma += Form("%.1f",ptTriggerMax); filenameSigma += "PtaFrom";
 624   filenameSigma += Form("%.1f",ptAssociatedMin); filenameSigma += "To";
 625   filenameSigma += Form("%.1f",ptAssociatedMax);  filenameSigma += ".txt";
 626   filenameSkewness += Form("%.1f",ptTriggerMin); filenameSkewness += "To";
 627   filenameSkewness += Form("%.1f",ptTriggerMax); filenameSkewness += "PtaFrom";
 628   filenameSkewness += Form("%.1f",ptAssociatedMin); filenameSkewness += "To";
 629   filenameSkewness += Form("%.1f",ptAssociatedMax);
 630   filenameSkewness += ".txt";
 631   filenameKurtosis += Form("%.1f",ptTriggerMin); filenameKurtosis += "To";
 632   filenameKurtosis += Form("%.1f",ptTriggerMax); filenameKurtosis += "PtaFrom";
 633   filenameKurtosis += Form("%.1f",ptAssociatedMin); filenameKurtosis += "To";
 634   filenameKurtosis += Form("%.1f",ptAssociatedMax);
 635   filenameKurtosis += ".txt";
 636
 637   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
 638   meanLatex = "#mu = ";
 639   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
 640   meanLatex += " #pm ";
 641   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
 642   ofstream fileMean(filenameMean.Data(),ios::app);
 643   fileMean << gCentrality << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
 644   fileMean.close();
 645
 646   rmsLatex = "#sigma = ";
 647   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
 648   rmsLatex += " #pm ";
 649   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
 650   ofstream fileSigma(filenameSigma.Data(),ios::app);
 651   fileSigma << gCentrality << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
 652   fileSigma.close();
 653
 654   skewnessLatex = "S = ";
 655   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
 656   skewnessLatex += " #pm ";
 657   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
 658   ofstream fileSkewness(filenameSkewness.Data(),ios::app);
 659   fileSkewness << gCentrality << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
 660   fileSkewness.close();
 661
 662   kurtosisLatex = "K = ";
 663   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
 664   kurtosisLatex += " #pm ";
 665   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
 666   ofstream fileKurtosis(filenameKurtosis.Data(),ios::app);
 667   fileKurtosis << gCentrality << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
 668   fileKurtosis.close();
 669
 670   Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
 671   Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
 672   Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
 673   Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
 674
 675   latexInfo1->DrawLatex(0.18,0.88,centralityLatex.Data());
 676   latexInfo1->DrawLatex(0.18,0.82,psiLatex.Data());
 677   latexInfo1->DrawLatex(0.18,0.76,pttLatex.Data());
 678   latexInfo1->DrawLatex(0.18,0.70,ptaLatex.Data());
 679
 680   TLatex *latexResults = new TLatex();
 681   latexResults->SetNDC();
 682   latexResults->SetTextSize(0.04);
 683   latexResults->SetTextColor(1);
 684   latexResults->DrawLatex(0.6,0.88,meanLatex.Data());
 685   latexResults->DrawLatex(0.6,0.82,rmsLatex.Data());
 686   latexResults->DrawLatex(0.6,0.76,skewnessLatex.Data());
 687   latexResults->DrawLatex(0.6,0.70,kurtosisLatex.Data());
 688
 689 }