PWGCF/EBYE/macros/drawCorrelationFunctionPsiChargeIndependent.C

   1 const Int_t numberOfCentralityBins = 13;
   2 TString centralityArray[numberOfCentralityBins] = {"0-80","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3","92-8500"};
   3
   4 const Int_t gRebin = 1;
   5
   6 void drawCorrelationFunctionPsiChargeIndependent(const char* filename = "AnalysisResults.root",
   7                                                  Int_t gCentrality = 1,
   8                                                  Int_t gBit = -1,
   9                                                  const char* gCentralityEstimator = 0x0,
  10                                                  Bool_t kShowShuffled = kFALSE,
  11                                                  Bool_t kShowMixed = kTRUE,
  12                                                  Double_t psiMin = -0.5,
  13                                                  Double_t psiMax = 3.5,
  14                                                  Double_t vertexZMin = -10.,
  15                                                  Double_t vertexZMax = 10.,
  16                                                  Double_t ptTriggerMin = -1.,
  17                                                  Double_t ptTriggerMax = -1.,
  18                                                  Double_t ptAssociatedMin = -1.,
  19                                                  Double_t ptAssociatedMax = -1.,
  20                                                  Bool_t normToTrig = kFALSE,
  21                                                  Int_t rebinEta = 1,
  22                                                  Int_t rebinPhi = 1) {
  23   //Macro that draws the correlation functions from the balance function
  24   //analysis vs the reaction plane
  25   //Author: Panos.Christakoglou@nikhef.nl
  26   gROOT->LoadMacro("~/SetPlotStyle.C");
  27   SetPlotStyle();
  28   gStyle->SetPalette(1,0);
  29
  30   //Load the PWG2ebye library
  31   gSystem->Load("libANALYSIS");
  32   gSystem->Load("libANALYSISalice");
  33   gSystem->Load("libEventMixing");
  34   gSystem->Load("libCORRFW");
  35   gSystem->Load("libPWGTools");
  36   gSystem->Load("libPWGCFebye");
  37
  38   //Prepare the objects and return them
  39   TList *list = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  40   TList *listShuffled = NULL;
  41   if(kShowShuffled) listShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  42   TList *listMixed = NULL;
  43   if(kShowMixed) listMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  44
  45   if(!list) {
  46     Printf("The TList object was not created");
  47     return;
  48   }
  49   else
  50     draw(list,listShuffled,listMixed,
  51          gCentralityEstimator,gCentrality,psiMin,psiMax,vertexZMin,vertexZMax,
  52          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig,rebinEta,rebinPhi);
  53 }
  54
  55 //______________________________________________________//
  56 TList *GetListOfObjects(const char* filename,
  57                         Int_t gCentrality,
  58                         Int_t gBit,
  59                         const char *gCentralityEstimator,
  60                         Int_t kData = 1) {
  61   //Get the TList objects (QA, bf, bf shuffled)
  62   TList *listBF = 0x0;
  63
  64   //Open the file
  65   TFile *f = TFile::Open(filename,"UPDATE");
  66   if((!f)||(!f->IsOpen())) {
  67     Printf("The file %s is not found. Aborting...",filename);
  68     return listBF;
  69   }
  70   //f->ls();
  71
  72   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  73   if(!dir) {
  74     Printf("The TDirectoryFile is not found. Aborting...",filename);
  75     return listBF;
  76   }
  77   //dir->ls();
  78
  79   TString listBFName;
  80   if(kData == 0) {
  81     //cout<<"no shuffling - no mixing"<<endl;
  82     listBFName = "listBFPsi_";
  83   }
  84   else if(kData == 1) {
  85     //cout<<"shuffling - no mixing"<<endl;
  86     listBFName = "listBFPsiShuffled_";
  87   }
  88   else if(kData == 2) {
  89     //cout<<"no shuffling - mixing"<<endl;
  90     listBFName = "listBFPsiMixed_";
  91   }
  92   listBFName += centralityArray[gCentrality-1];
  93   if(gBit > -1) {
  94     listBFName += "_Bit"; listBFName += gBit; }
  95   if(gCentralityEstimator) {
  96     listBFName += "_"; listBFName += gCentralityEstimator;}
  97
  98   // histograms were already retrieved (in first iteration)
  99   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
 100     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
 101   }
 102
 103   // histograms were not yet retrieved (this is the first iteration)
 104   else{
 105
 106     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 107     cout<<"======================================================="<<endl;
 108     cout<<"List name (control): "<<listBFName.Data()<<endl;
 109     cout<<"List name: "<<listBF->GetName()<<endl;
 110     //listBF->ls();
 111
 112     //Get the histograms
 113     TString histoName;
 114     if(kData == 0)
 115       histoName = "fHistP";
 116     else if(kData == 1)
 117       histoName = "fHistP_shuffle";
 118     else if(kData == 2)
 119       histoName = "fHistP";
 120     if(gCentralityEstimator)
 121       histoName += gCentralityEstimator;
 122     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 123     if(!fHistP) {
 124       Printf("fHistP %s not found!!!",histoName.Data());
 125       break;
 126     }
 127     fHistP->FillParent(); fHistP->DeleteContainers();
 128
 129     if(kData == 0)
 130       histoName = "fHistN";
 131     if(kData == 1)
 132       histoName = "fHistN_shuffle";
 133     if(kData == 2)
 134       histoName = "fHistN";
 135     if(gCentralityEstimator)
 136       histoName += gCentralityEstimator;
 137     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 138     if(!fHistN) {
 139       Printf("fHistN %s not found!!!",histoName.Data());
 140       break;
 141     }
 142     fHistN->FillParent(); fHistN->DeleteContainers();
 143
 144     if(kData == 0)
 145       histoName = "fHistPN";
 146     if(kData == 1)
 147       histoName = "fHistPN_shuffle";
 148     if(kData == 2)
 149       histoName = "fHistPN";
 150     if(gCentralityEstimator)
 151       histoName += gCentralityEstimator;
 152     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 153     if(!fHistPN) {
 154       Printf("fHistPN %s not found!!!",histoName.Data());
 155       break;
 156     }
 157     fHistPN->FillParent(); fHistPN->DeleteContainers();
 158
 159     if(kData == 0)
 160       histoName = "fHistNP";
 161     if(kData == 1)
 162       histoName = "fHistNP_shuffle";
 163     if(kData == 2)
 164       histoName = "fHistNP";
 165     if(gCentralityEstimator)
 166       histoName += gCentralityEstimator;
 167     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 168     if(!fHistNP) {
 169       Printf("fHistNP %s not found!!!",histoName.Data());
 170       break;
 171     }
 172     fHistNP->FillParent(); fHistNP->DeleteContainers();
 173
 174     if(kData == 0)
 175       histoName = "fHistPP";
 176     if(kData == 1)
 177       histoName = "fHistPP_shuffle";
 178     if(kData == 2)
 179       histoName = "fHistPP";
 180     if(gCentralityEstimator)
 181       histoName += gCentralityEstimator;
 182     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 183     if(!fHistPP) {
 184       Printf("fHistPP %s not found!!!",histoName.Data());
 185       break;
 186     }
 187     fHistPP->FillParent(); fHistPP->DeleteContainers();
 188
 189     if(kData == 0)
 190       histoName = "fHistNN";
 191     if(kData == 1)
 192       histoName = "fHistNN_shuffle";
 193     if(kData == 2)
 194       histoName = "fHistNN";
 195     if(gCentralityEstimator)
 196       histoName += gCentralityEstimator;
 197     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 198     if(!fHistNN) {
 199       Printf("fHistNN %s not found!!!",histoName.Data());
 200       break;
 201     }
 202     fHistNN->FillParent(); fHistNN->DeleteContainers();
 203
 204     dir->cd();
 205     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 206
 207   }// first iteration
 208
 209   f->Close();
 210
 211   return listBF;
 212 }
 213
 214 //______________________________________________________//
 215 void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
 216           const char *gCentralityEstimator,
 217           Int_t gCentrality, Double_t psiMin, Double_t psiMax,
 218           Double_t vertexZMin,
 219           Double_t vertexZMax,
 220           Double_t ptTriggerMin, Double_t ptTriggerMax,
 221           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 222           Bool_t normToTrig, Int_t rebinEta, Int_t rebinPhi) {
 223   //Draws the correlation functions for every centrality bin
 224   //(+-), (-+), (++), (--)
 225   AliTHn *hP = NULL;
 226   AliTHn *hN = NULL;
 227   AliTHn *hPN = NULL;
 228   AliTHn *hNP = NULL;
 229   AliTHn *hPP = NULL;
 230   AliTHn *hNN = NULL;
 231
 232   TString gHistPname = "fHistP";
 233   if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 234   TString gHistNname = "fHistN";
 235   if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 236   TString gHistPNname = "fHistPN";
 237   if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 238   TString gHistNPname = "fHistNP";
 239   if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 240   TString gHistPPname = "fHistPP";
 241   if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 242   TString gHistNNname = "fHistNN";
 243   if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 244
 245   hP = (AliTHn*) list->FindObject(gHistPname.Data());
 246   hN = (AliTHn*) list->FindObject(gHistNname.Data());
 247   hPN = (AliTHn*) list->FindObject(gHistPNname.Data());
 248   hNP = (AliTHn*) list->FindObject(gHistNPname.Data());
 249   hPP = (AliTHn*) list->FindObject(gHistPPname.Data());
 250   hNN = (AliTHn*) list->FindObject(gHistNNname.Data());
 251
 252   //Create the AliBalancePsi object and fill it with the AliTHn objects
 253   AliBalancePsi *b = new AliBalancePsi();
 254   b->SetHistNp(hP);
 255   b->SetHistNn(hN);
 256   b->SetHistNpn(hPN);
 257   b->SetHistNnp(hNP);
 258   b->SetHistNpp(hPP);
 259   b->SetHistNnn(hNN);
 260
 261   //balance function shuffling
 262   AliTHn *hPShuffled = NULL;
 263   AliTHn *hNShuffled = NULL;
 264   AliTHn *hPNShuffled = NULL;
 265   AliTHn *hNPShuffled = NULL;
 266   AliTHn *hPPShuffled = NULL;
 267   AliTHn *hNNShuffled = NULL;
 268   if(listBFShuffled) {
 269     //listBFShuffled->ls();
 270
 271     gHistPname = "fHistP_shuffle";
 272     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 273     gHistNname = "fHistN_shuffle";
 274     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 275     gHistPNname = "fHistPN_shuffle";
 276     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 277     gHistNPname = "fHistNP_shuffle";
 278     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 279     gHistPPname = "fHistPP_shuffle";
 280     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 281     gHistNNname = "fHistNN_shuffle";
 282     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 283
 284     hPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPname.Data());
 285     hPShuffled->SetName("gHistPShuffled");
 286     hNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNname.Data());
 287     hNShuffled->SetName("gHistNShuffled");
 288     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPNname.Data());
 289     hPNShuffled->SetName("gHistPNShuffled");
 290     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNPname.Data());
 291     hNPShuffled->SetName("gHistNPShuffled");
 292     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPPname.Data());
 293     hPPShuffled->SetName("gHistPPShuffled");
 294     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNNname.Data());
 295     hNNShuffled->SetName("gHistNNShuffled");
 296
 297     AliBalancePsi *bShuffled = new AliBalancePsi();
 298     bShuffled->SetHistNp(hPShuffled);
 299     bShuffled->SetHistNn(hNShuffled);
 300     bShuffled->SetHistNpn(hPNShuffled);
 301     bShuffled->SetHistNnp(hNPShuffled);
 302     bShuffled->SetHistNpp(hPPShuffled);
 303     bShuffled->SetHistNnn(hNNShuffled);
 304   }
 305
 306   //balance function mixing
 307   AliTHn *hPMixed = NULL;
 308   AliTHn *hNMixed = NULL;
 309   AliTHn *hPNMixed = NULL;
 310   AliTHn *hNPMixed = NULL;
 311   AliTHn *hPPMixed = NULL;
 312   AliTHn *hNNMixed = NULL;
 313
 314   if(listBFMixed) {
 315     //listBFMixed->ls();
 316
 317     gHistPname = "fHistP";
 318     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 319     gHistNname = "fHistN";
 320     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 321     gHistPNname = "fHistPN";
 322     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 323     gHistNPname = "fHistNP";
 324     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 325     gHistPPname = "fHistPP";
 326     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 327     gHistNNname = "fHistNN";
 328     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 329     hPMixed = (AliTHn*) listBFMixed->FindObject(gHistPname.Data());
 330     hPMixed->SetName("gHistPMixed");
 331     hNMixed = (AliTHn*) listBFMixed->FindObject(gHistNname.Data());
 332     hNMixed->SetName("gHistNMixed");
 333     hPNMixed = (AliTHn*) listBFMixed->FindObject(gHistPNname.Data());
 334     hPNMixed->SetName("gHistPNMixed");
 335     hNPMixed = (AliTHn*) listBFMixed->FindObject(gHistNPname.Data());
 336     hNPMixed->SetName("gHistNPMixed");
 337     hPPMixed = (AliTHn*) listBFMixed->FindObject(gHistPPname.Data());
 338     hPPMixed->SetName("gHistPPMixed");
 339     hNNMixed = (AliTHn*) listBFMixed->FindObject(gHistNNname.Data());
 340     hNNMixed->SetName("gHistNNMixed");
 341
 342     AliBalancePsi *bMixed = new AliBalancePsi();
 343     bMixed->SetHistNp(hPMixed);
 344     bMixed->SetHistNn(hNMixed);
 345     bMixed->SetHistNpn(hPNMixed);
 346     bMixed->SetHistNnp(hNPMixed);
 347     bMixed->SetHistNpp(hPPMixed);
 348     bMixed->SetHistNnn(hNNMixed);
 349   }
 350
 351   TH2D *gHist[6];
 352
 353   TCanvas *c[6];
 354   TString histoTitle, pngName;
 355
 356   // all charges together
 357   histoTitle = "(+-) | Centrality: ";
 358   histoTitle += centralityArray[gCentrality-1];
 359   histoTitle += "%";
 360   if((psiMin == -0.5)&&(psiMax == 0.5))
 361     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 362   else if((psiMin == 0.5)&&(psiMax == 1.5))
 363     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 364   else if((psiMin == 1.5)&&(psiMax == 2.5))
 365     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 366   else
 367     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 368
 369   gHist[0] = b->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 370   if(rebinEta > 1 || rebinPhi > 1){
 371     gHist[0]->Rebin2D(rebinEta,rebinPhi);
 372     gHist[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 373   }
 374   gHist[0]->GetYaxis()->SetTitleOffset(1.5);
 375   gHist[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 376   gHist[0]->SetTitle(histoTitle.Data());
 377   c[0] = new TCanvas("c0","",0,0,600,500);
 378   c[0]->SetFillColor(10); c[0]->SetHighLightColor(10);
 379   gHist[0]->DrawCopy("surf1fb");
 380   gPad->SetTheta(30); // default is 30
 381   //gPad->SetPhi(130); // default is 30
 382   gPad->SetPhi(-60); // default is 30
 383   gPad->Update();
 384   pngName = "DeltaPhiDeltaEta.Centrality";
 385   pngName += centralityArray[gCentrality-1];
 386   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 387   pngName += ".PositiveNegative.png";
 388   //c[0]->SaveAs(pngName.Data());
 389
 390   if(listBFShuffled) {
 391     histoTitle = "(+-) shuffled | Centrality: ";
 392     histoTitle += centralityArray[gCentrality-1];
 393     histoTitle += "%";
 394     if((psiMin == -0.5)&&(psiMax == 0.5))
 395       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 396     else if((psiMin == 0.5)&&(psiMax == 1.5))
 397       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 398     else if((psiMin == 1.5)&&(psiMax == 2.5))
 399       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 400     else
 401       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 402
 403     gHist[1] = bShuffled->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 404     if(rebinEta > 1 || rebinPhi > 1){
 405       gHist[1]->Rebin2D(rebinEta,rebinPhi);
 406       gHist[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 407     }
 408     gHist[1]->GetYaxis()->SetTitleOffset(1.5);
 409     gHist[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 410     gHist[1]->SetTitle(histoTitle.Data());
 411     c[1] = new TCanvas("c1","",0,100,600,500);
 412     c[1]->SetFillColor(10);
 413     c[1]->SetHighLightColor(10);
 414     gHist[1]->DrawCopy("surf1fb");
 415     gPad->SetTheta(30); // default is 30
 416     //gPad->SetPhi(130); // default is 30
 417     gPad->SetPhi(-60); // default is 30
 418     gPad->Update();
 419     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 420     pngName += centralityArray[gCentrality-1];
 421     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 422     pngName += ".PositiveNegative.png";
 423     //c[1]->SaveAs(pngName.Data());
 424   }
 425
 426   if(listBFMixed) {
 427     histoTitle = "(+-) mixed | Centrality: ";
 428     histoTitle += centralityArray[gCentrality-1];
 429     histoTitle += "%";
 430     if((psiMin == -0.5)&&(psiMax == 0.5))
 431       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 432     else if((psiMin == 0.5)&&(psiMax == 1.5))
 433       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 434     else if((psiMin == 1.5)&&(psiMax == 2.5))
 435       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 436     else
 437       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 438
 439     // if normalization to trigger then do not divide Event mixing by number of trigger particles
 440     gHist[2] = bMixed->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 441     if(rebinEta > 1 || rebinPhi > 1){
 442       gHist[2]->Rebin2D(rebinEta,rebinPhi);
 443       gHist[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 444     }
 445
 446     // normalization to 1 at (0,0) --> Jan Fietes method
 447     if(normToTrig){
 448       Double_t mixedNorm = gHist[2]->Integral(gHist[2]->GetXaxis()->FindBin(0-10e-5),gHist[2]->GetXaxis()->FindBin(0+10e-5),gHist[2]->GetNbinsY()/2 + 1,gHist[2]->GetNbinsY());
 449       mixedNorm /= 0.5 * gHist[2]->GetNbinsY()*(gHist[2]->GetXaxis()->FindBin(0.01) - gHist[2]->GetXaxis()->FindBin(-0.01) + 1);
 450
 451       // finite bin correction
 452       Double_t binWidthEta = gHist[2]->GetXaxis()->GetBinWidth(gHist[2]->GetNbinsX());
 453       Double_t maxEta      = gHist[2]->GetXaxis()->GetBinUpEdge(gHist[2]->GetNbinsX());
 454
 455       Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
 456       //Printf("Finite bin correction: %f", finiteBinCorrection);
 457       mixedNorm /= finiteBinCorrection;
 458
 459       gHist[2]->Scale(1./mixedNorm);
 460     }
 461
 462     gHist[2]->GetYaxis()->SetTitleOffset(1.5);
 463     gHist[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 464     gHist[2]->SetTitle(histoTitle.Data());
 465     c[2] = new TCanvas("c2","",0,200,600,500);
 466     c[2]->SetFillColor(10);
 467     c[2]->SetHighLightColor(10);
 468     gHist[2]->DrawCopy("surf1fb");
 469     gPad->SetTheta(30); // default is 30
 470     //gPad->SetPhi(130); // default is 30
 471     gPad->SetPhi(-60); // default is 30
 472     gPad->Update();
 473     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 474     pngName += centralityArray[gCentrality-1];
 475     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 476     pngName += ".PositiveNegative.png";
 477     //c[2]->SaveAs(pngName.Data());
 478
 479     //Correlation function (+-)
 480     gHist[3] = b->GetCorrelationFunction("ALL",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed,normToTrig);
 481     gHist[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 482     gHist[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 483     c[3] = new TCanvas("c3","",0,300,600,500);
 484     c[3]->SetFillColor(10);
 485     c[3]->SetHighLightColor(10);
 486     gHist[3]->DrawCopy("surf1fb");
 487     gPad->SetTheta(30); // default is 30
 488     //gPad->SetPhi(130); // default is 30
 489     gPad->SetPhi(-60); // default is 30
 490     gPad->Update();
 491     pngName = "CorrelationFunction.Centrality";
 492     pngName += centralityArray[gCentrality-1];
 493     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 494     pngName += ".PositiveNegative.png";
 495     //c[3]->SaveAs(pngName.Data());
 496
 497     //Correlation function subtracted awayside (+-)
 498     gHist[4] = dynamic_cast<TH2D *>(gHist[3]->Clone()); // this will be used to imitate twice the away-side
 499     gHist[4]->GetXaxis()->SetRangeUser(-1.5,1.5);
 500     gHist[4]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 501     gHist[5] = dynamic_cast<TH2D *>(gHist[3]->Clone()); // this will be the subtracted one
 502     gHist[5]->GetXaxis()->SetRangeUser(-1.5,1.5);
 503     gHist[5]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 504
 505     //prepare the double away side histo
 506     for(Int_t ix = 0; ix < gHist[4]->GetNbinsX(); ix++ ){
 507       for(Int_t iy = 0; iy < gHist[4]->GetNbinsX(); iy++ ){
 508         if(iy<gHist[4]->GetNbinsY()/2) gHist[4]->SetBinContent(ix+1,iy+1,gHist[4]->GetBinContent(ix+1,iy+1+gHist[4]->GetNbinsY()/2));
 509       }
 510     }
 511     gHist[5]->Add(gHist[4],-1);
 512
 513     c[4] = new TCanvas("c3","",0,300,600,500);
 514     c[4]->SetFillColor(10);
 515     c[4]->SetHighLightColor(10);
 516     gHist[5]->DrawCopy("surf1fb");
 517     gPad->SetTheta(30); // default is 30
 518     //gPad->SetPhi(130); // default is 30
 519     gPad->SetPhi(-60); // default is 30
 520     gPad->Update();
 521   }
 522
 523   //Write to output file
 524   TString newFileName = "correlationFunctionChargeIndependent.Centrality";
 525   newFileName += gCentrality; newFileName += ".Psi";
 526   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 527   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 528   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 529   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 530   else newFileName += "All.PttFrom";
 531   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 532   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 533   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 534   newFileName += Form("%.1f",ptAssociatedMax);
 535   newFileName += ".root";
 536   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 537   gHist[0]->SetName("gHistRaw"); gHist[0]->Write();
 538   if(listBFShuffled) {
 539     gHist[1]->SetName("gHistShuffled"); gHist[1]->Write();
 540   }
 541   if(listBFMixed) {
 542     gHist[2]->SetName("gHistMixed"); gHist[2]->Write();
 543
 544     gHist[3]->SetName("gHistCorrelationFunctions");gHist[3]->Write();
 545     gHist[4]->SetName("gHistCorrelationFunctionsAwaySide"); gHist[4]->Write();
 546     gHist[5]->SetName("gHistCorrelationFunctionsSubtracted"); gHist[5]->Write();
 547   }
 548   newFile->Close();
 549
 550   // // some cleaning
 551   // for(Int_t i = 0; i < 6; i++){
 552
 553   //   if(!listBFShuffled && i == 1) continue;
 554   //   if(!listBFMixed && (i == 2 || i == 3 || i == 4 || i == 5)) continue;
 555
 556   //   if(gHist[i]) delete gHist[i];
 557
 558   //   if(c[i]) delete c[i];
 559   // }
 560
 561   delete hP;
 562   delete hN;
 563   delete hPP;
 564   delete hPN;
 565   delete hNP;
 566   delete hNN;
 567
 568   delete hPMixed;
 569   delete hNMixed;
 570   delete hPPMixed;
 571   delete hPNMixed;
 572   delete hNPMixed;
 573   delete hNNMixed;
 574
 575   delete hPShuffled;
 576   delete hNShuffled;
 577   delete hPPShuffled;
 578   delete hPNShuffled;
 579   delete hNPShuffled;
 580   delete hNNShuffled;
 581
 582 }
 583
 584 //____________________________________________________________//
 585 void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
 586                               Int_t gTrainID = 208,
 587                               Int_t gCentrality = 1,
 588                               Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 589                               Double_t ptTriggerMin = -1.,
 590                               Double_t ptTriggerMax = -1.,
 591                               Double_t ptAssociatedMin = -1.,
 592                               Double_t ptAssociatedMax = -1.) {
 593   //Macro that draws the charge dependent correlation functions
 594   //for each centrality bin for the different pT of trigger and
 595   //associated particles
 596   //Author: Panos.Christakoglou@nikhef.nl
 597   TGaxis::SetMaxDigits(3);
 598
 599   //Get the input file
 600   TString filename = "PbPb/"; filename += lhcPeriod;
 601   filename +="/Train"; filename += gTrainID;
 602   filename +="/Centrality"; filename += gCentrality;
 603   filename += "/correlationFunction.Centrality";
 604   filename += gCentrality; filename += ".Psi";
 605   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 606   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 607   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 608   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 609   else filename += "All.PttFrom";
 610   filename += Form("%.1f",ptTriggerMin); filename += "To";
 611   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 612   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 613   filename += Form("%.1f",ptAssociatedMax);
 614   filename += ".root";
 615
 616   //Open the file
 617   TFile *f = TFile::Open(filename.Data());
 618   if((!f)||(!f->IsOpen())) {
 619     Printf("The file %s is not found. Aborting...",filename);
 620     return listBF;
 621   }
 622   //f->ls();
 623
 624   //Latex
 625   TString centralityLatex = "Centrality: ";
 626   centralityLatex += centralityArray[gCentrality-1];
 627   centralityLatex += "%";
 628
 629   TString psiLatex;
 630   if((psiMin == -0.5)&&(psiMax == 0.5))
 631     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 632   else if((psiMin == 0.5)&&(psiMax == 1.5))
 633     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 634   else if((psiMin == 1.5)&&(psiMax == 2.5))
 635     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 636   else
 637     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 638
 639   TString pttLatex = Form("%.1f",ptTriggerMin);
 640   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
 641   pttLatex += " GeV/c";
 642
 643   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 644   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 645   ptaLatex += " GeV/c";
 646
 647   TLatex *latexInfo1 = new TLatex();
 648   latexInfo1->SetNDC();
 649   latexInfo1->SetTextSize(0.045);
 650   latexInfo1->SetTextColor(1);
 651
 652   TString pngName;
 653
 654   //============================================================//
 655   //Get the +- correlation function
 656   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
 657   gHistPN->SetStats(kFALSE);
 658   gHistPN->SetTitle("");
 659   gHistPN->GetXaxis()->SetRangeUser(-1.45,1.45);
 660   gHistPN->GetXaxis()->CenterTitle();
 661   gHistPN->GetXaxis()->SetTitleOffset(1.2);
 662   gHistPN->GetYaxis()->CenterTitle();
 663   gHistPN->GetYaxis()->SetTitleOffset(1.2);
 664   gHistPN->GetZaxis()->SetTitleOffset(1.2);
 665   TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
 666   cPN->SetFillColor(10); cPN->SetHighLightColor(10);
 667   cPN->SetLeftMargin(0.15);
 668   gHistPN->DrawCopy("surf1fb");
 669   gPad->SetTheta(30); // default is 30
 670   gPad->SetPhi(-60); // default is 30
 671   gPad->Update();
 672
 673   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
 674   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
 675   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
 676   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
 677
 678   pngName = "CorrelationFunctionChargeIndependent.Centrality";
 679   pngName += centralityArray[gCentrality-1];
 680   pngName += ".Psi";
 681   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 682   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 683   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 684   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
 685   else pngName += "All.PttFrom";
 686   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 687   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 688   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 689   pngName += Form("%.1f",ptAssociatedMax);
 690   pngName += ".PositiveNegative.png";
 691   cPN->SaveAs(pngName.Data());
 692   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
 693                           ptTriggerMin,ptTriggerMax,
 694                           ptAssociatedMin, ptAssociatedMax,gHistPN);
 695   //============================================================//
 696   //Get the -+ correlation function
 697   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
 698   gHistNP->SetStats(kFALSE);
 699   gHistNP->SetTitle("");
 700   gHistNP->GetXaxis()->SetRangeUser(-1.45,1.45);
 701   gHistNP->GetXaxis()->CenterTitle();
 702   gHistNP->GetXaxis()->SetTitleOffset(1.2);
 703   gHistNP->GetYaxis()->CenterTitle();
 704   gHistNP->GetYaxis()->SetTitleOffset(1.2);
 705   gHistNP->GetZaxis()->SetTitleOffset(1.2);
 706   TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
 707   cNP->SetFillColor(10); cNP->SetHighLightColor(10);
 708   cNP->SetLeftMargin(0.15);
 709   gHistNP->DrawCopy("surf1fb");
 710   gPad->SetTheta(30); // default is 30
 711   gPad->SetPhi(-60); // default is 30
 712   gPad->Update();
 713
 714   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
 715   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
 716   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
 717   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
 718
 719   pngName = "CorrelationFunction.Centrality";
 720   pngName += centralityArray[gCentrality-1];
 721   pngName += ".Psi";
 722   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 723   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 724   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 725   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
 726   else pngName += "All.PttFrom";
 727   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 728   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 729   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 730   pngName += Form("%.1f",ptAssociatedMax);
 731   pngName += ".NegativePositive.png";
 732   cNP->SaveAs(pngName.Data());
 733
 734   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
 735                           ptTriggerMin,ptTriggerMax,
 736                           ptAssociatedMin, ptAssociatedMax,gHistNP);
 737   //============================================================//
 738   //Get the ++ correlation function
 739   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
 740   gHistPP->SetStats(kFALSE);
 741   gHistPP->SetTitle("");
 742   gHistPP->GetXaxis()->SetRangeUser(-1.45,1.45);
 743   gHistPP->GetXaxis()->CenterTitle();
 744   gHistPP->GetXaxis()->SetTitleOffset(1.2);
 745   gHistPP->GetYaxis()->CenterTitle();
 746   gHistPP->GetYaxis()->SetTitleOffset(1.2);
 747   gHistPP->GetZaxis()->SetTitleOffset(1.2);
 748   TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
 749   cPP->SetFillColor(10); cPP->SetHighLightColor(10);
 750   cPP->SetLeftMargin(0.15);
 751   gHistPP->DrawCopy("surf1fb");
 752   gPad->SetTheta(30); // default is 30
 753   gPad->SetPhi(-60); // default is 30
 754   gPad->Update();
 755
 756   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
 757   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
 758   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
 759   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
 760
 761   pngName = "CorrelationFunction.Centrality";
 762   pngName += centralityArray[gCentrality-1];
 763   pngName += ".Psi";
 764   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 765   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 766   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 767   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
 768   else pngName += "All.PttFrom";
 769   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 770   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 771   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 772   pngName += Form("%.1f",ptAssociatedMax);
 773   pngName += ".PositivePositive.png";
 774   cPP->SaveAs(pngName.Data());
 775
 776   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
 777                           ptTriggerMin,ptTriggerMax,
 778                           ptAssociatedMin, ptAssociatedMax,gHistPP);
 779   //============================================================//
 780   //Get the -- correlation function
 781   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
 782   gHistNN->SetStats(kFALSE);
 783   gHistNN->SetTitle("");
 784   gHistNN->GetXaxis()->SetRangeUser(-1.45,1.45);
 785   gHistNN->GetXaxis()->CenterTitle();
 786   gHistNN->GetXaxis()->SetTitleOffset(1.2);
 787   gHistNN->GetYaxis()->CenterTitle();
 788   gHistNN->GetYaxis()->SetTitleOffset(1.2);
 789   gHistNN->GetZaxis()->SetTitleOffset(1.2);
 790   TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
 791   cNN->SetFillColor(10); cNN->SetHighLightColor(10);
 792   cNN->SetLeftMargin(0.15);
 793   gHistNN->DrawCopy("surf1fb");
 794   gPad->SetTheta(30); // default is 30
 795   gPad->SetPhi(-60); // default is 30
 796   gPad->Update();
 797
 798   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
 799   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
 800   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
 801   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
 802
 803   pngName = "CorrelationFunction.Centrality";
 804   pngName += centralityArray[gCentrality-1];
 805   pngName += ".Psi";
 806   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 807   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 808   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 809   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
 810   else pngName += "All.PttFrom";
 811   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 812   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 813   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 814   pngName += Form("%.1f",ptAssociatedMax);
 815   pngName += ".NegativeNegative.png";
 816   cNN->SaveAs(pngName.Data());
 817
 818   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
 819                           ptTriggerMin,ptTriggerMax,
 820                           ptAssociatedMin, ptAssociatedMax,gHistNN);
 821 }
 822
 823 // //____________________________________________________________//
 824 // void fitCorrelationFunctions(Int_t gCentrality = 1,
 825 //                           Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 826 //                           Double_t ptTriggerMin = -1.,
 827 //                           Double_t ptTriggerMax = -1.,
 828 //                           Double_t ptAssociatedMin = -1.,
 829 //                           Double_t ptAssociatedMax = -1.,
 830 //                           TH2D *gHist) {
 831
 832 //   cout<<"FITTING FUNCTION (MW style)"<<endl;
 833
 834 //   //near side peak(s): [1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) +
 835 //   //                        TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4])))
 836 //   //away side peak(s): [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power((y/[9]),2)),[10])) +
 837 //   //                        TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power((y/[9]),2)),[10])))
 838 //   //flow contribution (v1 up to v4): 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))
 839
 840
 841
 842 //   TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4]))) + [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10]))) + 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
 843 //   gFitFunction->SetName("gFitFunction");
 844
 845
 846 //   //Normalization
 847 //   gFitFunction->SetParName(0,"N1");
 848 //   //near side peak(s)
 849 //   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
 850 //   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(2,0.0);
 851 //   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
 852 //   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
 853 //   gFitFunction->SetParName(5,"Offset_{near side}"); gFitFunction->FixParameter(5,0.0);
 854 //   gFitFunction->SetParName(6,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(6,0.0);
 855
 856 //   //away side peak(s)
 857 //   gFitFunction->SetParName(7,"N_{near side}");gFitFunction->FixParameter(7,0.0);
 858 //   gFitFunction->SetParName(8,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(8,0.0);
 859 //   gFitFunction->SetParName(9,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(9,0.0);
 860 //   gFitFunction->SetParName(10,"Exponent_{near side}"); gFitFunction->FixParameter(10,1.0);
 861 //   gFitFunction->SetParName(11,"Offset_{near side}"); gFitFunction->FixParameter(11,0.0);
 862 //   gFitFunction->SetParName(12,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(12,0.0);
 863
 864 //   //flow contribution
 865 //   gFitFunction->SetParName(13,"N_{flow}"); gFitFunction->SetParameter(13,0.2);
 866 //   gFitFunction->SetParName(14,"V1"); gFitFunction->SetParameter(14,0.005);
 867 //   gFitFunction->SetParName(15,"V2"); gFitFunction->SetParameter(15,0.1);
 868 //   gFitFunction->SetParName(16,"V3"); gFitFunction->SetParameter(16,0.05);
 869 //   gFitFunction->SetParName(17,"V4"); gFitFunction->SetParameter(17,0.005);
 870
 871 //   // flow parameters
 872 //   Double_t fNV = 0.;
 873 //   Double_t fV1 = 0.;
 874 //   Double_t fV2 = 0.;
 875 //   Double_t fV3 = 0.;
 876 //   Double_t fV4 = 0.;
 877
 878 //   //Fitting the correlation function (first the edges to extract flow)
 879 //   gHist->Fit("gFitFunction","nm","",1.0,1.6);
 880
 881 //   fNV += gFitFunction->GetParameter(13);
 882 //   fV1 += gFitFunction->GetParameter(14);
 883 //   fV2 += gFitFunction->GetParameter(15);
 884 //   fV3 += gFitFunction->GetParameter(16);
 885 //   fV4 += gFitFunction->GetParameter(17);
 886
 887 //   gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
 888
 889 //   fNV += gFitFunction->GetParameter(13);
 890 //   fV1 += gFitFunction->GetParameter(14);
 891 //   fV2 += gFitFunction->GetParameter(15);
 892 //   fV3 += gFitFunction->GetParameter(16);
 893 //   fV4 += gFitFunction->GetParameter(17);
 894
 895 //   // Now fit the whole with fixed flow
 896 //   gFitFunction->FixParameter(13,fNV/2.);
 897 //   gFitFunction->FixParameter(14,fV1/2.);
 898 //   gFitFunction->FixParameter(15,fV2/2.);
 899 //   gFitFunction->FixParameter(16,fV3/2.);
 900 //   gFitFunction->FixParameter(17,fV4/2.);
 901
 902 //   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
 903 //   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
 904 //   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
 905 //   gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
 906 //   gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
 907
 908 //   gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(1,0.3);
 909 //   gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
 910 //   gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
 911 //   gFitFunction->ReleaseParameter(11);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
 912 //   gFitFunction->ReleaseParameter(12);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
 913
 914 //   gHist->Fit("gFitFunction","nm");
 915
 916
 917 //   //Cloning the histogram
 918 //   TString histoName = gHist->GetName(); histoName += "Fit";
 919 //   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());
 920 //   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
 921 //   gHistResidual->SetName("gHistResidual");
 922 //   gHistResidual->Sumw2();
 923
 924 //   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
 925 //     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
 926 //       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
 927 //     }
 928 //   }
 929 //   gHistResidual->Add(gHistFit,-1);
 930
 931 //   //Write to output file
 932 //   TString newFileName = "correlationFunctionFit";
 933 //   if(histoName.Contains("PN")) newFileName += "PN";
 934 //   else if(histoName.Contains("NP")) newFileName += "NP";
 935 //   else if(histoName.Contains("PP")) newFileName += "PP";
 936 //   else if(histoName.Contains("NN")) newFileName += "NN";
 937 //   newFileName += ".Centrality";
 938 //   newFileName += gCentrality; newFileName += ".Psi";
 939 //   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 940 //   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 941 //   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 942 //   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 943 //   else newFileName += "All.PttFrom";
 944 //   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 945 //   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 946 //   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 947 //   newFileName += Form("%.1f",ptAssociatedMax);
 948 //   newFileName += ".root";
 949 //   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 950 //   gHist->Write();
 951 //   gHistFit->Write();
 952 //   gHistResidual->Write();
 953 //   gFitFunction->Write();
 954 //   newFile->Close();
 955
 956
 957 // }
 958
 959 //____________________________________________________________//
 960 void fitCorrelationFunctions(Int_t gCentrality = 1,
 961                              Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 962                              Double_t ptTriggerMin = -1.,
 963                              Double_t ptTriggerMax = -1.,
 964                              Double_t ptAssociatedMin = -1.,
 965                              Double_t ptAssociatedMax = -1.,
 966                              TH2D *gHist) {
 967
 968   cout<<"FITTING FUNCTION"<<endl;
 969
 970   //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
 971   //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
 972   //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
 973   //wing structures: [11]*TMath::Power(x,2)
 974   //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))
 975   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.);
 976   gFitFunction->SetName("gFitFunction");
 977
 978
 979   //Normalization
 980   gFitFunction->SetParName(0,"N1");
 981   //near side peak
 982   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
 983   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
 984   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
 985   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
 986   //away side ridge
 987   gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
 988   gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
 989   gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
 990   //longitudinal ridge
 991   gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
 992   gFitFunction->FixParameter(8,0.0);
 993   gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
 994   gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
 995   //wing structures
 996   gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
 997   //flow contribution
 998   gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);gFitFunction->SetParLimits(12,0.0,10);
 999   gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);gFitFunction->SetParLimits(13,0.0,10);
1000   gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);gFitFunction->SetParLimits(14,0.0,10);
1001   gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);gFitFunction->SetParLimits(15,0.0,10);
1002   gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);gFitFunction->SetParLimits(16,0.0,10);
1003   gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
1004
1005   // flow parameters
1006   Double_t fNV = 0.;
1007   Double_t fV1 = 0.;
1008   Double_t fV2 = 0.;
1009   Double_t fV3 = 0.;
1010   Double_t fV4 = 0.;
1011
1012   //Fitting the correlation function (first the edges to extract flow)
1013   gHist->Fit("gFitFunction","nm","",1.0,1.6);
1014
1015   fNV += gFitFunction->GetParameter(12);
1016   fV1 += gFitFunction->GetParameter(13);
1017   fV2 += gFitFunction->GetParameter(14);
1018   fV3 += gFitFunction->GetParameter(15);
1019   fV4 += gFitFunction->GetParameter(16);
1020
1021   gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
1022
1023   fNV += gFitFunction->GetParameter(12);
1024   fV1 += gFitFunction->GetParameter(13);
1025   fV2 += gFitFunction->GetParameter(14);
1026   fV3 += gFitFunction->GetParameter(15);
1027   fV4 += gFitFunction->GetParameter(16);
1028
1029   // Now fit the whole with fixed flow
1030   gFitFunction->FixParameter(12,fNV/2.);
1031   gFitFunction->FixParameter(13,fV1/2.);
1032   gFitFunction->FixParameter(14,fV2/2.);
1033   gFitFunction->FixParameter(15,fV3/2.);
1034   gFitFunction->FixParameter(16,fV4/2.);
1035
1036   gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
1037   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
1038   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
1039   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
1040   //gFitFunction->ReleaseParameter(4);gFitFunction->SetParameter(4,1.0);gFitFunction->SetParLimits(4,0.0,2.0);
1041   gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,1.0);//gFitFunction->SetParLimits(5,0.0,10);
1042   gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.5);//gFitFunction->SetParLimits(6,0.0,10);
1043   //gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(7,1.0);gFitFunction->SetParLimits(7,0.0,2.0);
1044   gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(8,0.05);
1045   gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(9,0.6);gFitFunction->SetParLimits(9,0.1,10.0);
1046   //gFitFunction->ReleaseParameter(10);gFitFunction->SetParameter(10,1.0);gFitFunction->SetParLimits(10,0.0,2.0);
1047   gFitFunction->ReleaseParameter(17);gFitFunction->SetParameter(17,0.7);gFitFunction->SetParLimits(17,0.0,0.9);
1048
1049   gHist->Fit("gFitFunction","nm");
1050
1051
1052   //Cloning the histogram
1053   TString histoName = gHist->GetName(); histoName += "Fit";
1054   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());
1055   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1056   gHistResidual->SetName("gHistResidual");
1057   gHistResidual->Sumw2();
1058
1059   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1060     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1061       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1062     }
1063   }
1064   gHistResidual->Add(gHistFit,-1);
1065
1066   //Write to output file
1067   TString newFileName = "correlationFunctionFit";
1068   if(histoName.Contains("PN")) newFileName += "PN";
1069   else if(histoName.Contains("NP")) newFileName += "NP";
1070   else if(histoName.Contains("PP")) newFileName += "PP";
1071   else if(histoName.Contains("NN")) newFileName += "NN";
1072   newFileName += ".Centrality";
1073   newFileName += gCentrality; newFileName += ".Psi";
1074   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1075   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1076   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1077   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1078   else newFileName += "All.PttFrom";
1079   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1080   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1081   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1082   newFileName += Form("%.1f",ptAssociatedMax);
1083   newFileName += ".root";
1084   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1085   gHist->Write();
1086   gHistFit->Write();
1087   gHistResidual->Write();
1088   gFitFunction->Write();
1089   newFile->Close();
1090
1091
1092 }
1093
1094
1095
1096 // //____________________________________________________________//
1097 // void fitCorrelationFunctions(Int_t gCentrality = 1,
1098 //                           Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1099 //                           Double_t ptTriggerMin = -1.,
1100 //                           Double_t ptTriggerMax = -1.,
1101 //                           Double_t ptAssociatedMin = -1.,
1102 //                           Double_t ptAssociatedMax = -1.,
1103 //                           TH2D *gHist) {
1104
1105 //   cout<<"FITTING FUNCTION (HOUSTON)"<<endl;
1106
1107 //   // Fit Function
1108 //    //x axis = delta_eta
1109 //    //y axis = delta_phi
1110 //                                                                                                                                                                                         //  [9]*exp(-1*pow(((x/[10])^2 + (y/[10])^2),0.5)) Hyper expo
1111 //   TF2 *fit1 = new TF2("fit1","[0] + [1]*cos(y) + [2]*cos(2*y) + [3]*cos(3*y) + [4]*cos(4*y) +[5]*cos(5*y)+ [6]*exp(-0.5*pow(((x/[7])^2 + (y/[8])^2),[11])) + [6]*exp(-0.5*pow(((x/[7])^2 + ((y-6.283)/[8])^2),[11]))+ [9]*exp(-1*((x/[10])^2 + (y/[10])^2)) ",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
1112
1113
1114
1115 //   Double_t Parameters[] = {0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.3,0.3,1.0,0.1,0.1};
1116
1117 //   fit1->SetParameters(Parameters);  // input pars from macro arguments
1118
1119 //   fit1->SetParName(0,"offset");
1120 //   fit1->SetParName(1,"v1");
1121 //   fit1->SetParName(2,"v2");
1122 //   fit1->SetParName(3,"v3");
1123 //   fit1->SetParName(4,"v4");
1124 //   fit1->SetParName(5,"v5");
1125 //   fit1->SetParName(6,"Ridgeamp");
1126 //   fit1->SetParName(7,"Ridgesigx");
1127 //   fit1->SetParName(8,"Ridgesigy");
1128 //   fit1->SetParName(9,"Expoamp");
1129 //   fit1->SetParName(10,"Exposig");
1130 //   fit1->SetParName(11,"Gausspara");
1131
1132
1133 //   //Fit Parameter Ranges
1134 //   fit1->SetParLimits(0,-2.0,2.0);   //offset
1135 //   fit1->SetParLimits(1,-1.0,0.1);   //v1
1136 //   fit1->SetParLimits(2,-1.6,0.9);   //v2
1137 //   fit1->SetParLimits(3,0.0,0.5);    //v3
1138 //   fit1->SetParLimits(4,0.0,0.9);    //v4
1139 //   fit1->SetParLimits(5,0.0,0.9);    //v5
1140 //   fit1->SetParLimits(6,0.0,1.5);    //Ridgeamp
1141 //   fit1->SetParLimits(7,0.1,3.0);    //Ridgesigx
1142 //   fit1->SetParLimits(8,0.1,2.0);    //Ridgesigy
1143 //   fit1->SetParLimits(9,0.0,6.0);      //Expoamp
1144 //   fit1->SetParLimits(10,0.05,0.5); //Exposig
1145 //   fit1->SetParLimits(11,0.0,2.0);    //Gausspara
1146
1147 //   //Fitting the correlation function
1148 //   gHist->Fit("fit1","nm");
1149
1150 //   //Cloning the histogram
1151 //   TString histoName = gHist->GetName(); histoName += "Fit";
1152 //   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());
1153 //   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1154 //   gHistResidual->SetName("gHistResidual");
1155 //   gHistResidual->Sumw2();
1156
1157 //   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1158 //     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1159 //       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,fit1->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1160 //     }
1161 //   }
1162 //   gHistResidual->Add(gHistFit,-1);
1163
1164 //   //Write to output file
1165 //   TString newFileName = "correlationFunctionFit";
1166 //   if(histoName.Contains("PN")) newFileName += "PN";
1167 //   else if(histoName.Contains("NP")) newFileName += "NP";
1168 //   else if(histoName.Contains("PP")) newFileName += "PP";
1169 //   else if(histoName.Contains("NN")) newFileName += "NN";
1170 //   newFileName += ".Centrality";
1171 //   newFileName += gCentrality; newFileName += ".Psi";
1172 //   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1173 //   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1174 //   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1175 //   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1176 //   else newFileName += "All.PttFrom";
1177 //   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1178 //   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1179 //   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1180 //   newFileName += Form("%.1f",ptAssociatedMax);
1181 //   newFileName += ".root";
1182 //   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1183 //   gHist->Write();
1184 //   gHistFit->Write();
1185 //   gHistResidual->Write();
1186 //   fit1->Write();
1187 //   newFile->Close();
1188
1189
1190 // }