PWGCF/EBYE/macros/drawCorrelationFunctionPsiChargeIndependent.C

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