PWGCF/EBYE/macros/drawBalanceFunction2DPsi.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 const Int_t gRebin = 1;
   5 void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
   6                               Int_t gCentrality = 1,
   7                               Int_t gBit = -1,
   8                               const char* gCentralityEstimator = 0x0,
   9                               Bool_t kShowShuffled = kFALSE,
  10                               Bool_t kShowMixed = kTRUE,
  11                               Double_t psiMin = -0.5, Double_t psiMax = 0.5,
  12                               Double_t ptTriggerMin = -1.,
  13                               Double_t ptTriggerMax = -1.,
  14                               Double_t ptAssociatedMin = -1.,
  15                               Double_t ptAssociatedMax = -1.,
  16                               Bool_t k2pMethod = kFALSE) {
  17   //Macro that draws the BF distributions for each centrality bin
  18   //for reaction plane dependent analysis
  19   //Author: Panos.Christakoglou@nikhef.nl
  20   //Load the PWG2ebye library
  21   gSystem->Load("libANALYSIS.so");
  22   gSystem->Load("libANALYSISalice.so");
  23   gSystem->Load("libEventMixing.so");
  24   gSystem->Load("libCORRFW.so");
  25   gSystem->Load("libPWGTools.so");
  26   gSystem->Load("libPWGCFebye.so");
  27
  28   gROOT->LoadMacro("~/SetPlotStyle.C");
  29   SetPlotStyle();
  30   gStyle->SetPalette(1,0);
  31
  32   //Prepare the objects and return them
  33   TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  34   TList *listBFShuffled = NULL;
  35   if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  36   TList *listBFMixed = NULL;
  37   if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  38   if(!listBF) {
  39     Printf("The TList object was not created");
  40     return;
  41   }
  42   else
  43     draw(listBF,listBFShuffled,listBFMixed,gCentrality,psiMin,psiMax,
  44          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
  45          k2pMethod);
  46 }
  47
  48 //______________________________________________________//
  49 TList *GetListOfObjects(const char* filename,
  50                         Int_t gCentrality,
  51                         Int_t gBit,
  52                         const char *gCentralityEstimator,
  53                         Int_t kData = 1) {
  54   //Get the TList objects (QA, bf, bf shuffled)
  55   TList *listBF = 0x0;
  56
  57   //Open the file
  58   TFile *f = TFile::Open(filename,"UPDATE");
  59   if((!f)||(!f->IsOpen())) {
  60     Printf("The file %s is not found. Aborting...",filename);
  61     return listBF;
  62   }
  63   //f->ls();
  64
  65   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  66   if(!dir) {
  67     Printf("The TDirectoryFile is not found. Aborting...",filename);
  68     return listBF;
  69   }
  70   //dir->ls();
  71
  72   TString listBFName;
  73   if(kData == 0) {
  74     //cout<<"no shuffling - no mixing"<<endl;
  75     listBFName = "listBFPsi_";
  76   }
  77   else if(kData == 1) {
  78     //cout<<"shuffling - no mixing"<<endl;
  79     listBFName = "listBFPsiShuffled_";
  80   }
  81   else if(kData == 2) {
  82     //cout<<"no shuffling - mixing"<<endl;
  83     listBFName = "listBFPsiMixed_";
  84   }
  85   listBFName += centralityArray[gCentrality-1];
  86   if(gBit > -1) {
  87     listBFName += "_Bit"; listBFName += gBit; }
  88   if(gCentralityEstimator) {
  89     listBFName += "_"; listBFName += gCentralityEstimator;}
  90
  91   // histograms were already retrieved (in first iteration)
  92   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
  93     //cout<<"second iteration"<<endl;
  94     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
  95   }
  96
  97   // histograms were not yet retrieved (this is the first iteration)
  98   else{
  99
 100     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 101     cout<<"======================================================="<<endl;
 102     cout<<"List name (check): "<<listBFName.Data()<<endl;
 103     cout<<"List name: "<<listBF->GetName()<<endl;
 104     //listBF->ls();
 105
 106     //Get the histograms
 107     TString histoName;
 108     if(kData == 0)
 109       histoName = "fHistPV0M";
 110     else if(kData == 1)
 111       histoName = "fHistP_shuffleV0M";
 112     else if(kData == 2)
 113       histoName = "fHistPV0M";
 114     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 115     if(!fHistP) {
 116       Printf("fHistP %s not found!!!",histoName.Data());
 117       break;
 118     }
 119     fHistP->FillParent(); fHistP->DeleteContainers();
 120
 121     if(kData == 0)
 122       histoName = "fHistNV0M";
 123     if(kData == 1)
 124       histoName = "fHistN_shuffleV0M";
 125     if(kData == 2)
 126       histoName = "fHistNV0M";
 127     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 128     if(!fHistN) {
 129       Printf("fHistN %s not found!!!",histoName.Data());
 130       break;
 131     }
 132     fHistN->FillParent(); fHistN->DeleteContainers();
 133
 134     if(kData == 0)
 135       histoName = "fHistPNV0M";
 136     if(kData == 1)
 137       histoName = "fHistPN_shuffleV0M";
 138     if(kData == 2)
 139       histoName = "fHistPNV0M";
 140     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 141     if(!fHistPN) {
 142       Printf("fHistPN %s not found!!!",histoName.Data());
 143       break;
 144     }
 145     fHistPN->FillParent(); fHistPN->DeleteContainers();
 146
 147     if(kData == 0)
 148       histoName = "fHistNPV0M";
 149     if(kData == 1)
 150       histoName = "fHistNP_shuffleV0M";
 151     if(kData == 2)
 152       histoName = "fHistNPV0M";
 153     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 154     if(!fHistNP) {
 155       Printf("fHistNP %s not found!!!",histoName.Data());
 156       break;
 157     }
 158     fHistNP->FillParent(); fHistNP->DeleteContainers();
 159
 160     if(kData == 0)
 161       histoName = "fHistPPV0M";
 162     if(kData == 1)
 163       histoName = "fHistPP_shuffleV0M";
 164     if(kData == 2)
 165       histoName = "fHistPPV0M";
 166     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 167     if(!fHistPP) {
 168       Printf("fHistPP %s not found!!!",histoName.Data());
 169       break;
 170     }
 171     fHistPP->FillParent(); fHistPP->DeleteContainers();
 172
 173     if(kData == 0)
 174       histoName = "fHistNNV0M";
 175     if(kData == 1)
 176       histoName = "fHistNN_shuffleV0M";
 177     if(kData == 2)
 178       histoName = "fHistNNV0M";
 179     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 180     if(!fHistNN) {
 181       Printf("fHistNN %s not found!!!",histoName.Data());
 182       break;
 183     }
 184     fHistNN->FillParent(); fHistNN->DeleteContainers();
 185
 186     dir->cd();
 187     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 188
 189   }// first iteration
 190
 191   f->Close();
 192
 193   return listBF;
 194 }
 195
 196 //______________________________________________________//
 197 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 198           Int_t gCentrality, Double_t psiMin, Double_t psiMax,
 199           Double_t ptTriggerMin, Double_t ptTriggerMax,
 200           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 201           Bool_t k2pMethod = kFALSE) {
 202   //balance function
 203   AliTHn *hP = NULL;
 204   AliTHn *hN = NULL;
 205   AliTHn *hPN = NULL;
 206   AliTHn *hNP = NULL;
 207   AliTHn *hPP = NULL;
 208   AliTHn *hNN = NULL;
 209   //listBF->ls();
 210   //Printf("=================");
 211   hP = (AliTHn*) listBF->FindObject("fHistPV0M");
 212   hP->SetName("gHistP");
 213   hN = (AliTHn*) listBF->FindObject("fHistNV0M");
 214   hN->SetName("gHistN");
 215   hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
 216   hPN->SetName("gHistPN");
 217   hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
 218   hNP->SetName("gHistNP");
 219   hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
 220   hPP->SetName("gHistPP");
 221   hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
 222   hNN->SetName("gHistNN");
 223
 224   AliBalancePsi *b = new AliBalancePsi();
 225   b->SetHistNp(hP);
 226   b->SetHistNn(hN);
 227   b->SetHistNpn(hPN);
 228   b->SetHistNnp(hNP);
 229   b->SetHistNpp(hPP);
 230   b->SetHistNnn(hNN);
 231
 232   //balance function shuffling
 233   AliTHn *hPShuffled = NULL;
 234   AliTHn *hNShuffled = NULL;
 235   AliTHn *hPNShuffled = NULL;
 236   AliTHn *hNPShuffled = NULL;
 237   AliTHn *hPPShuffled = NULL;
 238   AliTHn *hNNShuffled = NULL;
 239   if(listBFShuffled) {
 240     //listBFShuffled->ls();
 241
 242     hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
 243     hPShuffled->SetName("gHistPShuffled");
 244     hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
 245     hNShuffled->SetName("gHistNShuffled");
 246     hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
 247     hPNShuffled->SetName("gHistPNShuffled");
 248     hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
 249     hNPShuffled->SetName("gHistNPShuffled");
 250     hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
 251     hPPShuffled->SetName("gHistPPShuffled");
 252     hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
 253     hNNShuffled->SetName("gHistNNShuffled");
 254
 255     AliBalancePsi *bShuffled = new AliBalancePsi();
 256     bShuffled->SetHistNp(hPShuffled);
 257     bShuffled->SetHistNn(hNShuffled);
 258     bShuffled->SetHistNpn(hPNShuffled);
 259     bShuffled->SetHistNnp(hNPShuffled);
 260     bShuffled->SetHistNpp(hPPShuffled);
 261     bShuffled->SetHistNnn(hNNShuffled);
 262   }
 263
 264   //balance function mixing
 265   AliTHn *hPMixed = NULL;
 266   AliTHn *hNMixed = NULL;
 267   AliTHn *hPNMixed = NULL;
 268   AliTHn *hNPMixed = NULL;
 269   AliTHn *hPPMixed = NULL;
 270   AliTHn *hNNMixed = NULL;
 271
 272   if(listBFMixed) {
 273     //listBFMixed->ls();
 274
 275     hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
 276     hPMixed->SetName("gHistPMixed");
 277     hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
 278     hNMixed->SetName("gHistNMixed");
 279     hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
 280     hPNMixed->SetName("gHistPNMixed");
 281     hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
 282     hNPMixed->SetName("gHistNPMixed");
 283     hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
 284     hPPMixed->SetName("gHistPPMixed");
 285     hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
 286     hNNMixed->SetName("gHistNNMixed");
 287
 288     AliBalancePsi *bMixed = new AliBalancePsi();
 289     bMixed->SetHistNp(hPMixed);
 290     bMixed->SetHistNn(hNMixed);
 291     bMixed->SetHistNpn(hPNMixed);
 292     bMixed->SetHistNnp(hNPMixed);
 293     bMixed->SetHistNpp(hPPMixed);
 294     bMixed->SetHistNnn(hNNMixed);
 295   }
 296
 297   TH2D *gHistBalanceFunction;
 298   TH2D *gHistBalanceFunctionSubtracted;
 299   TH2D *gHistBalanceFunctionShuffled;
 300   TH2D *gHistBalanceFunctionMixed;
 301   TString histoTitle, pngName;
 302
 303   histoTitle = "Centrality: ";
 304   histoTitle += centralityArray[gCentrality-1];
 305   histoTitle += "%";
 306   if((psiMin == -0.5)&&(psiMax == 0.5))
 307     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 308   else if((psiMin == 0.5)&&(psiMax == 1.5))
 309     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 310   else if((psiMin == 1.5)&&(psiMax == 2.5))
 311     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 312   else
 313     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 314
 315   if(k2pMethod)
 316     if(bMixed)
 317       gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 318     else{
 319       cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 320       return;
 321     }
 322   else
 323     gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 324   gHistBalanceFunction->SetTitle(histoTitle.Data());
 325   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 326   gHistBalanceFunction->SetName("gHistBalanceFunction");
 327
 328   if(listBFShuffled) {
 329
 330     if(k2pMethod)
 331       if(bMixed)
 332         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 333       else{
 334         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 335         return;
 336       }
 337     else
 338       gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 339     gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
 340     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 341     gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 342   }
 343
 344   if(listBFMixed) {
 345     if(k2pMethod)
 346       if(bMixed)
 347         gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 348       else{
 349         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 350         return;
 351       }
 352     else
 353       gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 354     gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
 355     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 356     gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 357
 358     gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
 359     gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 360     gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 361     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 362     gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 363   }
 364
 365   //Draw the results
 366   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 367   c1->SetFillColor(10);
 368   c1->SetHighLightColor(10);
 369   c1->SetLeftMargin(0.15);
 370   gHistBalanceFunction->DrawCopy("lego2");
 371   gPad->SetTheta(30); // default is 30
 372   gPad->SetPhi(-60); // default is 30
 373   gPad->Update();
 374   TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
 375   c1a->SetFillColor(10);
 376   c1a->SetHighLightColor(10);
 377   c1a->SetLeftMargin(0.15);
 378   gHistBalanceFunction->DrawCopy("colz");
 379
 380   if(listBFShuffled) {
 381     TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
 382     c2->SetFillColor(10);
 383     c2->SetHighLightColor(10);
 384     c2->SetLeftMargin(0.15);
 385     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 386     gPad->SetTheta(30); // default is 30
 387     gPad->SetPhi(-60); // default is 30
 388     gPad->Update();
 389     TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
 390     c2a->SetFillColor(10);
 391     c2a->SetHighLightColor(10);
 392     c2a->SetLeftMargin(0.15);
 393     gHistBalanceFunctionShuffled->DrawCopy("colz");
 394   }
 395
 396   if(listBFMixed) {
 397     TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
 398     c3->SetFillColor(10);
 399     c3->SetHighLightColor(10);
 400     c3->SetLeftMargin(0.15);
 401     gHistBalanceFunctionMixed->DrawCopy("lego2");
 402     gPad->SetTheta(30); // default is 30
 403     gPad->SetPhi(-60); // default is 30
 404     gPad->Update();
 405     TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
 406     c3a->SetFillColor(10);
 407     c3a->SetHighLightColor(10);
 408     c3a->SetLeftMargin(0.15);
 409     gHistBalanceFunctionMixed->DrawCopy("colz");
 410
 411     TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
 412     c4->SetFillColor(10);
 413     c4->SetHighLightColor(10);
 414     c4->SetLeftMargin(0.15);
 415     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 416     gPad->SetTheta(30); // default is 30
 417     gPad->SetPhi(-60); // default is 30
 418     gPad->Update();
 419     TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
 420     c4a->SetFillColor(10);
 421     c4a->SetHighLightColor(10);
 422     c4a->SetLeftMargin(0.15);
 423     gHistBalanceFunctionSubtracted->DrawCopy("colz");
 424
 425     fitbalanceFunction(gCentrality, psiMin = -0.5, psiMax,
 426                        ptTriggerMin, ptTriggerMax,
 427                        ptAssociatedMin, ptAssociatedMax,
 428                        gHistBalanceFunctionSubtracted);
 429   }
 430
 431   TString newFileName = "balanceFunction2D.Centrality";
 432   newFileName += gCentrality; newFileName += ".Psi";
 433   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 434   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 435   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 436   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 437   else newFileName += "All.PttFrom";
 438   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 439   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 440   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 441   newFileName += Form("%.1f",ptAssociatedMax);
 442   if(k2pMethod) newFileName += "_2pMethod";
 443   newFileName += ".root";
 444
 445   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 446   fOutput->cd();
 447   /*hP->Write(); hN->Write();
 448   hPN->Write(); hNP->Write();
 449   hPP->Write(); hNN->Write();
 450   hPShuffled->Write(); hNShuffled->Write();
 451   hPNShuffled->Write(); hNPShuffled->Write();
 452   hPPShuffled->Write(); hNNShuffled->Write();
 453   hPMixed->Write(); hNMixed->Write();
 454   hPNMixed->Write(); hNPMixed->Write();
 455   hPPMixed->Write(); hNNMixed->Write();*/
 456   gHistBalanceFunction->Write();
 457   if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
 458   if(listBFMixed) {
 459     gHistBalanceFunctionMixed->Write();
 460     gHistBalanceFunctionSubtracted->Write();
 461   }
 462   fOutput->Close();
 463 }
 464
 465 //____________________________________________________________//
 466 void fitbalanceFunction(Int_t gCentrality = 1,
 467                         Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 468                         Double_t ptTriggerMin = -1.,
 469                         Double_t ptTriggerMax = -1.,
 470                         Double_t ptAssociatedMin = -1.,
 471                         Double_t ptAssociatedMax = -1.,
 472                         TH2D *gHist) {
 473
 474   cout<<"FITTING FUNCTION"<<endl;
 475
 476   //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
 477   //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
 478   //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
 479   //wing structures: [11]*TMath::Power(x,2)
 480   //flow contribution (v1 up to v4): 2.*([12]*TMath::Cos(y) + [13]*TMath::Cos(2.*y) + [14]*TMath::Cos(3.*y) + [15]*TMath::Cos(4.*y))
 481   TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))+[5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[17])/[9]),2)),[10]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[17])/[9]),2)),[10]))+[11]*TMath::Power(x,2)+2.*[12]*([13]*TMath::Cos(y) + [14]*TMath::Cos(2.*y) + [15]*TMath::Cos(3.*y) + [16]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
 482   gFitFunction->SetName("gFitFunction");
 483
 484
 485   //Normalization
 486   gFitFunction->SetParName(0,"N1");
 487   //near side peak
 488   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
 489   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
 490   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
 491   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
 492   //away side ridge
 493   gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
 494   gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
 495   gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
 496   //longitudinal ridge
 497   gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
 498   gFitFunction->FixParameter(8,0.0);
 499   gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
 500   gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
 501   //wing structures
 502   gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
 503   //flow contribution
 504   gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->FixParameter(12,0.0);
 505   gFitFunction->SetParName(13,"V1"); gFitFunction->FixParameter(13,0.0);
 506   gFitFunction->SetParName(14,"V2"); gFitFunction->FixParameter(14,0.0);
 507   gFitFunction->SetParName(15,"V3"); gFitFunction->FixParameter(15,0.0);
 508   gFitFunction->SetParName(16,"V4"); gFitFunction->FixParameter(16,0.0);
 509   gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
 510
 511   // just release the near side peak
 512   gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
 513   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);gFitFunction->SetParLimits(1,0.0,100);
 514   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
 515   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
 516
 517   gHist->Fit("gFitFunction","nm");
 518
 519
 520   //Cloning the histogram
 521   TString histoName = gHist->GetName(); histoName += "Fit";
 522   TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
 523   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
 524   gHistResidual->SetName("gHistResidual");
 525   gHistResidual->Sumw2();
 526
 527   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
 528     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
 529       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
 530     }
 531   }
 532   gHistResidual->Add(gHistFit,-1);
 533
 534   //Write to output file
 535   TString newFileName = "balanceFunctionFit";
 536   newFileName += ".Centrality";
 537   newFileName += gCentrality; newFileName += ".Psi";
 538   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 539   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 540   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 541   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 542   else newFileName += "All.PttFrom";
 543   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 544   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 545   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 546   newFileName += Form("%.1f",ptAssociatedMax);
 547   newFileName += ".root";
 548   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 549   gHist->Write();
 550   gHistFit->Write();
 551   gHistResidual->Write();
 552   gFitFunction->Write();
 553   newFile->Close();
 554
 555
 556 }
 557
 558
 559
 560 //____________________________________________________________//
 561 void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
 562                  Int_t gCentrality = 1,
 563                  Bool_t kShowShuffled = kFALSE,
 564                  Bool_t kShowMixed = kFALSE,
 565                  Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 566                  Double_t ptTriggerMin = -1.,
 567                  Double_t ptTriggerMax = -1.,
 568                  Double_t ptAssociatedMin = -1.,
 569                  Double_t ptAssociatedMax = -1.) {
 570   //Macro that draws the BF distributions for each centrality bin
 571   //for reaction plane dependent analysis
 572   //Author: Panos.Christakoglou@nikhef.nl
 573   TGaxis::SetMaxDigits(3);
 574
 575   //Get the input file
 576   TString filename = lhcPeriod; filename +="/PttFrom";
 577   filename += Form("%.1f",ptTriggerMin); filename += "To";
 578   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 579   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 580   filename += Form("%.1f",ptAssociatedMax);
 581   filename += "/balanceFunction2D.Centrality";
 582   filename += gCentrality; filename += ".Psi";
 583   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 584   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 585   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 586   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 587   else filename += "All.PttFrom";
 588   filename += Form("%.1f",ptTriggerMin); filename += "To";
 589   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 590   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 591   filename += Form("%.1f",ptAssociatedMax);  filename += ".root";
 592
 593   //Open the file
 594   TFile *f = TFile::Open(filename.Data());
 595   if((!f)||(!f->IsOpen())) {
 596     Printf("The file %s is not found. Aborting...",filename);
 597     return listBF;
 598   }
 599   //f->ls();
 600
 601   //Raw balance function
 602   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 603   gHistBalanceFunction->SetStats(kFALSE);
 604   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 605   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 606   gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
 607   gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
 608   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 609   gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
 610   gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
 611   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
 612   gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 613
 614   //Shuffled balance function
 615   if(kShowShuffled) {
 616     TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 617     gHistBalanceFunctionShuffled->SetStats(kFALSE);
 618     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 619     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 620     gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
 621     gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
 622     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 623     gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
 624     gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
 625     gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
 626     gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 627   }
 628
 629   //Mixed balance function
 630   if(kShowMixed) {
 631     TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 632     gHistBalanceFunctionMixed->SetStats(kFALSE);
 633     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 634     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 635     gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
 636     gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
 637     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 638     gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
 639     gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
 640     gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
 641     gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 642   }
 643
 644   //Subtracted balance function
 645   if(kShowMixed) {
 646     TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 647     gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 648     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 649     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 650     gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
 651     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
 652     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 653     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
 654     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
 655     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
 656     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 657   }
 658
 659   TString pngName;
 660
 661   TString centralityLatex = "Centrality: ";
 662   centralityLatex += centralityArray[gCentrality-1];
 663   centralityLatex += "%";
 664
 665   TString psiLatex;
 666   if((psiMin == -0.5)&&(psiMax == 0.5))
 667     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 668   else if((psiMin == 0.5)&&(psiMax == 1.5))
 669     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 670   else if((psiMin == 1.5)&&(psiMax == 2.5))
 671     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 672   else
 673     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 674
 675   TString pttLatex = Form("%.1f",ptTriggerMin);
 676   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 677   pttLatex += " GeV/c";
 678
 679   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 680   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 681   ptaLatex += " GeV/c";
 682
 683   TLatex *latexInfo1 = new TLatex();
 684   latexInfo1->SetNDC();
 685   latexInfo1->SetTextSize(0.045);
 686   latexInfo1->SetTextColor(1);
 687
 688   //Draw the results
 689   TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
 690   c1->SetFillColor(10); c1->SetHighLightColor(10);
 691   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 692   gHistBalanceFunction->SetTitle("Raw balance function");
 693   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 694   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 695   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 696   gHistBalanceFunction->DrawCopy("lego2");
 697   gPad->SetTheta(30); // default is 30
 698   gPad->SetPhi(-60); // default is 30
 699   gPad->Update();
 700
 701   latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 702   latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 703   latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 704   latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 705
 706   if(kShowShuffled) {
 707     TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
 708     c2->SetFillColor(10); c2->SetHighLightColor(10);
 709     c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 710     gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
 711     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
 712     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 713     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 714     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 715     gPad->SetTheta(30); // default is 30
 716     gPad->SetPhi(-60); // default is 30
 717     gPad->Update();
 718
 719     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 720     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 721     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 722     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 723   }
 724
 725   if(kShowMixed) {
 726     TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
 727     c3->SetFillColor(10); c3->SetHighLightColor(10);
 728     c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
 729     gHistBalanceFunctionMixed->SetTitle("Mixed events");
 730     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
 731     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 732     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 733     gHistBalanceFunctionMixed->DrawCopy("lego2");
 734     gPad->SetTheta(30); // default is 30
 735     gPad->SetPhi(-60); // default is 30
 736     gPad->Update();
 737
 738     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 739     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 740     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 741     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 742   }
 743
 744   if(kShowMixed) {
 745     TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
 746     c4->SetFillColor(10); c4->SetHighLightColor(10);
 747     c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
 748     gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
 749     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 750     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 751     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 752     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 753     gPad->SetTheta(30); // default is 30
 754     gPad->SetPhi(-60); // default is 30
 755     gPad->Update();
 756
 757     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 758     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 759     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 760     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 761   }
 762 }