marcro updates
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Dec 2013 09:44:20 +0000 (09:44 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Dec 2013 09:44:20 +0000 (09:44 +0000)
PWGCF/Correlations/macros/dphicorrelations/correct.C
PWGCF/Correlations/macros/dphicorrelations/fit.C
PWGCF/Correlations/macros/dphicorrelations/pA.C
PWGCF/Correlations/macros/dphicorrelations/pA_PID.C
PWGCF/Correlations/macros/dphicorrelations/pA_PID_Fwd.C [new file with mode: 0644]

index 4862188..1973649 100644 (file)
@@ -493,6 +493,10 @@ void CompareEventHist(const char* fileName1, const char* fileName2, Int_t id, In
 
   TH1* hist1 = h->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
   TH1* hist2 = h2->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
+  
+//   hist1 = hist1->Rebin(hist2->GetNbinsX(), "new", hist2->GetXaxis()->GetXbins()->GetArray());
+  
+//   hist2 = hist2->Rebin(hist1->GetNbinsX(), "new", hist1->GetXaxis()->GetXbins()->GetArray());
  
   DrawRatio(hist1, hist2, "compare");
 }
@@ -852,6 +856,62 @@ void DrawExample(const char* fileName, const char* fileNamePbPbMix = 0)
   DrawALICELogo(0.75, 0.75, 0.9, 0.9);
 }
 
+void DrawExample2(const char* fileName)
+{
+  loadlibs();
+  
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+  hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
+  
+  gpTMin = 0.51;
+  gpTMax = 0.99;
+  SetupRanges(h);
+  SetupRanges(hMixed);
+
+  TH1* hist1 = 0;
+  GetSumOfRatios(h, hMixed, &hist1,  8, 0,  5, 0.51, 0.99);
+  
+  new TCanvas("c", "c", 800, 800);
+  gPad->SetLeftMargin(0.15);
+  hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  hist1->GetXaxis()->SetTitleOffset(1.5);
+  hist1->GetYaxis()->SetTitleOffset(2);
+  hist1->SetStats(kFALSE);
+  hist1->DrawCopy("SURF1"); 
+}  
+
+void CompareZVtxWeighting(const char* fileName, Int_t step = 8)
+{
+  gpTMin = 1.01;
+  gpTMax = 1.99;
+  
+  loadlibs();
+  
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+  hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
+  
+  SetupRanges(h);
+  SetupRanges(hMixed);
+  
+  TH1* hist1 = 0;
+  GetDistAndFlow(h, hMixed, &hist1,  0, step, 0,  10, 5.01, 9.99, 1, kTRUE, 0, kTRUE); 
+  hist1->Scale(1.0 / 0.972222);
+  
+  TH1* hist2 = 0;
+  GetSumOfRatios(h, hMixed, &hist2,  step, 0,  10, 5.01, 9.99, kTRUE);
+  
+  c = new TCanvas("c", "c", 1200, 400);
+  c->Divide(3, 1);
+  
+  c->cd(1); hist1->DrawCopy("SURF1");
+  c->cd(2); hist2->DrawCopy("SURF1");
+  
+  hist1->SetStats(0);
+  hist1->Divide(hist2);
+  
+  c->cd(3); hist1->DrawCopy("COLZ");
+}
+
 void DrawSameMixed(const char* fileName, const char* fileNamePbPbMix = 0, Int_t step = 6)
 {
   if (!fileNamePbPbMix)
@@ -904,7 +964,7 @@ void DrawSameMixed(const char* fileName, const char* fileNamePbPbMix = 0, Int_t
   hist2 = hist1;
   
   GetDistAndFlow(hMixed, 0, &hist1,  0, step, 0,  80, 2.01, 3.99, 1, kTRUE, 0, kTRUE); 
-
+  
 //   ((TH2*) hist1)->Rebin2D(2, 2);
   NormalizeToBinWidth(hist1);
   
@@ -2237,7 +2297,7 @@ void ModifyComposition(const char* fileNameData, const char* fileNameCorrections
     Printf("Largest deviation in region %d is %f", i, largestDeviation[i]);
 }    
 
-void MergeList(const char* prefix = "")
+void MergeList(const char* prefix = "", Bool_t copy = kFALSE)
 {
   loadlibs();
   gSystem->Load("libPWGflowBase");
@@ -2246,7 +2306,7 @@ void MergeList(const char* prefix = "")
   ifstream in;
   in.open("list");
 
-  TFileMerger m(0);
+  TFileMerger m(copy);
 
   TString line;
   while (in.good())
@@ -6641,7 +6701,7 @@ void* cacheIds[10];
 TH2* cacheMixed[10];
 
 Int_t gHistCount = 0;
-void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t* v2, Int_t step, Int_t centralityBegin, Int_t centralityEnd, Float_t ptBegin, Float_t ptEnd, Int_t twoD = 0, Bool_t equivMixedBin = kFALSE, Float_t* vn = 0, Bool_t scaleToPairs = kTRUE, Int_t stepMixed = -1)
+void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t* v2, Int_t step, Int_t centralityBegin, Int_t centralityEnd, Float_t ptBegin, Float_t ptEnd, Int_t twoD = 0, Bool_t equivMixedBin = kFALSE, Float_t* vn = 0, Bool_t scaleToPairs = kTRUE, Int_t stepMixed = -1, Bool_t useCentralityBinsDirectly = kFALSE)
 {
   h = (AliUEHistograms*) hVoid;
   hMixed = (AliUEHistograms*) hMixedVoid;
@@ -6652,11 +6712,16 @@ void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t* v2, Int_
   Int_t centralityBeginBin = 0;
   Int_t centralityEndBin = -1;
   
-  if (centralityEnd >= centralityBegin)
+  if (!useCentralityBinsDirectly && centralityEnd >= centralityBegin)
   {
     centralityBeginBin = h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(0.01 + centralityBegin);
     centralityEndBin = h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(-0.01 + centralityEnd);
   }
+  else if (useCentralityBinsDirectly)
+  {
+    centralityBeginBin = centralityBegin;
+    centralityEndBin = centralityEnd;
+  }
   
   // 2d same and mixed event
   TH2* sameTwoD  = h->GetUEHist(2)->GetUEHist(step, 0, ptBegin, ptEnd, centralityBeginBin, centralityEndBin, 1, kFALSE);
@@ -6810,7 +6875,9 @@ void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t* v2, Int_
   str2.Form("%.2f < p_{T,assoc} < %.2f", gpTMin - 0.01, gpTMax + 0.01);
     
   TString newTitle;
-  newTitle.Form("%s - %s - %d-%d%%", str.Data(), str2.Data(), centralityBegin, centralityEnd);
+  newTitle.Form("%s - %s - %d-%d", str.Data(), str2.Data(), centralityBegin, centralityEnd);
+  if (!useCentralityBinsDirectly)
+    newTitle += "%";
   (*hist)->SetTitle(newTitle);
   
   if (0 && hMixed)
@@ -7659,15 +7726,15 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
     
   Bool_t symmetrizePt = kFALSE;
   Int_t maxLeadingPt = 4;
-  Int_t maxAssocPt = 7;
-  if (0)
+  Int_t maxAssocPt = 6;
+  if (1)
   {
     //PbPb, NS peak shapes
-//     Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
-    Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 8.0, 15.0, 20.0 };
-    Float_t assocPtArr[] =     { 0.15, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0 };
+    Float_t leadingPtArr[] = { 1.0, 2.0, 3.0, 4.0, 8.0, 15.0, 20.0 };
+//     Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 8.0, 15.0, 20.0 };
+    Float_t assocPtArr[] =     { 0.15, 0.5, 1.0, 2.0, 3.0, 4.0, 8.0, 10.0, 12.0 };
   }
-  else if (1)
+  else if (0)
   {
     //Example for the Hadrons_Example wagon
     maxLeadingPt = 3;
@@ -7680,7 +7747,7 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
     //pA, trigger from all pT
     maxLeadingPt = 1;
     maxAssocPt = 10;
-    Float_t leadingPtArr[] =   { 0.5, 4.0 };
+    Float_t leadingPtArr[] =   { 0.3, 4.0 };
     Float_t assocPtArr[] =     { 0.15, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 0.5, 4.0 };
 //     symmetrizePt = kTRUE;
   }
@@ -7883,6 +7950,7 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
   }
   
   for (Int_t i=0; i<maxLeadingPt; i++)
+  {
     for (Int_t j=1; j<maxAssocPt; j++)
     {
       if(0){
@@ -7903,7 +7971,7 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
       SetupRanges(hMixed3);
 //       SetupRanges(hMixed3);
 
-      if(0)if (assocPtArr[j] >= leadingPtArr[i+leadingPtOffset])
+      if(1) if (assocPtArr[j] >= leadingPtArr[i+leadingPtOffset])
        continue;
   
       TH1* hist1 = 0;
@@ -7920,33 +7988,22 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
       
       Int_t histType = 1;
 
-      if (0)
+      if (1)
       {
        // PbPb
        Int_t step = 8;
+       Bool_t normalizePerTrigger = kFALSE; // don't do if histograms are to be merged -> Use MergeDPhiFiles below
       
-       GetSumOfRatios(h, hMixed, &hist1,     step, 0,  10, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-//     new TCanvas; hist1->DrawClone("SURF1");
-//     return;
+       GetSumOfRatios(h, hMixed, &hist1,  step, 0,   10, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+       GetSumOfRatios(h, hMixed, &hist5,  step, 10,  20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+       GetSumOfRatios(h, hMixed, &hist4,  step, 20,  40, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+       GetSumOfRatios(h, hMixed, &hist6,  step, 40,  60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+       GetSumOfRatios(h, hMixed, &hist2,  step, 60,  80, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+//     GetSumOfRatios(h, hMixed, &hist7,  step, 70,  80, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
+//     GetSumOfRatios(h, hMixed, &hist8,  step, 80,  90, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, normalizePerTrigger); 
 
-/*     GetDistAndFlow(h, hMixed, &hist2,  0, step, 0,   10,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs); 
-       new TCanvas; hist2->DrawClone("SURF1");
-       
-       hist1->Divide(hist2);
-       new TCanvas; hist1->DrawClone("SURF1");
-       
-       return;*/
-       
-//     new TCanvas; hist1->Draw("SURF1"); return;
-
-       GetSumOfRatios(h, hMixed, &hist5,  step, 10,  20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-       GetSumOfRatios(h, hMixed, &hist4,  step, 20,  40, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-       GetSumOfRatios(h, hMixed, &hist6,  step, 40,  60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-       GetSumOfRatios(h, hMixed, &hist2,  step, 60,  70, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-       step = 8;
        if (h2)
          GetSumOfRatios(h2, hMixed2, &hist3,  step, 0,  -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
-//     new TCanvas; hist3->Draw("SURF1"); return;
       }
       else if (0)
       {
@@ -7967,7 +8024,24 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
        if (h3)
          GetSumOfRatios(h3, hMixed3, &hist6,  step, 0,  -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
       }      
-      else if (1)
+      else if (0)
+      {
+       // pA, V0 as trigger particle detector
+       Int_t step = 6;
+      
+       h->GetUEHist(2)->SetSkipScaleMixedEvent(kTRUE);
+       GetSumOfRatios(h, hMixed, &hist1,  step,  0, 20, 1.01, 1.99, kTRUE); 
+       GetSumOfRatios(h, hMixed, &hist2,  step, 20, 40, 1.01, 1.99, kTRUE); 
+       GetSumOfRatios(h, hMixed, &hist4,  step, 40, 60, 1.01, 1.99, kTRUE); 
+       GetSumOfRatios(h, hMixed, &hist5,  step, 60, 100, 1.01, 1.99, kTRUE); 
+
+       if (h2)
+         GetSumOfRatios(h2, hMixed2, &hist3,  step, 0,  -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
+       
+       if (h3)
+         GetSumOfRatios(h3, hMixed3, &hist6,  step, 0,  -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
+      }      
+      else if (0)
       {
        // pp, MB
        Int_t step = 8;      
@@ -8025,14 +8099,14 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
        GetSumOfRatios(h, hMixed, &hist7,  10,  40, 60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
        GetSumOfRatios(h, hMixed, &hist8,  10,  60, 100, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
       }        
-      else if (1)
+      else if (0)
       {
        // pp, MB
        Int_t step = 8;
        
        GetSumOfRatios(h, hMixed, &hist1,  step,  0, 100, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
       }      
-      else if (1)
+      else if (0)
       {
        // pA, MC, validation binning, without vertex axis
        GetDistAndFlow(h, hMixed, &hist1,  0, 0,  0, 20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs); 
@@ -8046,6 +8120,23 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
       }        
       else if (1)
       {
+       // PbPb, MC, impact parameter binning
+       Int_t step = 0;
+       
+       Printf(">>>>>>>> Not using GetSumOfRatios!!!");
+
+       GetDistAndFlow(h, hMixed, &hist1,  0, step, 1,   6,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs, -1, kTRUE); 
+       GetDistAndFlow(h, hMixed, &hist5,  0, step, 7,  7, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs, -1, kTRUE); 
+       GetDistAndFlow(h, hMixed, &hist4,  0, step, 8,  9, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs, -1, kTRUE); 
+       GetDistAndFlow(h, hMixed, &hist6,  0, step, 10,  11, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs, -1, kTRUE); 
+       GetDistAndFlow(h, hMixed, &hist2,  0, step, 12,  13, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs, -1, kTRUE);
+
+       if (h2)
+         GetDistAndFlow(h2, hMixed2, &hist3,  0, step, 0, -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
+      }
+      else if (1)
+      {
+       // PbPb, MC
        Int_t step = 0;
        
        Printf(">>>>>>>> Not using GetSumOfRatios!!!");
@@ -8058,7 +8149,8 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
        GetDistAndFlow(h, hMixed, &hist6,  0, step, 40,  60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs); 
        GetDistAndFlow(h, hMixed, &hist2,  0, step, 60,  80, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
 //     step = 6;
-       GetDistAndFlow(h2, hMixed2, &hist3,  0, step, 0, -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
+       if (h2)
+         GetDistAndFlow(h2, hMixed2, &hist3,  0, step, 0, -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
 //     new TCanvas; hist3->DrawClone("SURF1");
 //     GetDistAndFlow(hMixed2, 0, &hist3,  0, step, 0, -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
 //     new TCanvas; hist3->DrawClone("SURF1");
@@ -8162,10 +8254,97 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix =
 //       return;
     }
     
+    TH1* triggers = h->GetUEHist(2)->GetTriggersAsFunctionOfMultiplicity(step, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01);
+    triggers->SetName(Form("triggers_%d", i));
+    TString str;
+    str.Form("%.1f < p_{T,trig} < %.1f", leadingPtArr[i], leadingPtArr[i+leadingPtOffset]);
+    triggers->SetTitle(str);
+
+    file = TFile::Open(outputFile, "UPDATE");
+    triggers->Write();
+    file->Close();    
+  }
+    
   delete h;
   delete hMixed;
 }
 
+void MergeDPhiFiles(const char* fileName, const char* fileName2, const char* target)
+{
+  // merges the dphi histograms (except the pp histogram at index "2") as well as the triggers
+  // then divides the dphi histograms by the respective number of triggers
+  file = TFile::Open(fileName);
+  file2 = TFile::Open(fileName2);
+
+  fileTarget = TFile::Open(target, "RECREATE");
+  fileTarget->Close();
+  
+  Int_t maxLeadingPt = 10;
+  Int_t maxAssocPt = 11;
+
+  Int_t nHists = 8;
+  for (Int_t i=0; i<maxLeadingPt; i++)
+  {
+    TH1* triggers = (TH1*) file->Get(Form("triggers_%d", i));
+    if (!triggers)
+      continue;
+    
+    TH1* triggers2 = (TH1*) ((file2) ? file2->Get(Form("triggers_%d", i)) : 0);
+    if (!triggers2)
+      Printf("WARNING: trigger %d missing", i);
+    
+    for (Int_t j=0; j<maxAssocPt; j++)
+    {
+      for (Int_t histId = 0; histId < nHists; histId++)
+      {
+       TH2* hist = (TH2*) file->Get(Form("dphi_%d_%d_%d", i, j, histId));
+       if (!hist)
+       {
+         TH2* hist2 = (TH2*) ((file2) ? file2->Get(Form("dphi_%d_%d_%d", i, j, histId)) : 0);
+         if (hist2)
+           Printf("WARNING: %d %d %d exists only in file2, not copied!");
+         continue;
+       }
+       
+       if (histId != 2) // don't merge pp
+       {
+         TString title(hist->GetTitle());
+         title.ReplaceAll("%", "");
+         tokens = title.Tokenize("-");
+         
+         Float_t centralityBegin = ((TObjString*) tokens->At(2))->String().Atoi();
+         Float_t centralityEnd = ((TObjString*) tokens->At(3))->String().Atoi();
+         
+         Double_t nTriggers = triggers->Integral(triggers->FindBin(centralityBegin + 0.001), triggers->FindBin(centralityEnd - 0.001));
+         Double_t nTriggers2 = 0;
+
+         TH2* hist2 = (TH2*) ((file2) ? file2->Get(Form("dphi_%d_%d_%d", i, j, histId)) : 0);
+         if (hist2 && triggers2) 
+         {
+           if (histId != 1 && histId != 5) // don't merge 60-80% and 40-60%
+           {
+             hist->Add(hist2);
+             nTriggers2 = triggers2->Integral(triggers2->FindBin(centralityBegin + 0.001), triggers2->FindBin(centralityEnd - 0.001));
+           }
+         }
+         else
+           Printf("WARNING: %d %d %d missing", i, j, histId);
+           
+         if (nTriggers + nTriggers2 > 0)
+           hist->Scale(1.0 / (nTriggers + nTriggers2));
+
+         Printf("%s %f %f %f %f", hist->GetTitle(), centralityBegin, centralityEnd, nTriggers, nTriggers2);
+       }
+
+       fileTarget = TFile::Open(target, "UPDATE");
+       hist->Write();
+       fileTarget->Close();
+      }
+    }
+  }
+}
+
 void ExtractMiniJetHistograms(const char* fileNamePbPb, Bool_t useMixed = kTRUE, const char* outputFile = "dphi_corr.root")
 {
   loadlibs();
@@ -8869,7 +9048,7 @@ void AnalyzeDeltaPhiEtaGap2D(const char* fileName, Int_t method)
 
 void RemoveWing(const char* fileName, const char* outputFile)
 {
-  // remove wing by flattening using the ratio of a flat line to the corr fct at phi = pi/2 +- 0.25 as fct of delta eta
+  // remove wing by flattening using the ratio of a flat line to the corr fct at phi = pi +- 1.5 as fct of delta eta
 
   file = TFile::Open(fileName);
   file2 = TFile::Open(outputFile, "RECREATE");
@@ -12167,25 +12346,32 @@ void PlotQA(const char* fileName, const char* tag = "")
   title.Form("QA_%d_%s", runNumber, tmp.Data());
   c->SetTitle(title);
   
-  gpTMin = 1.0;
+  gpTMin = 1.01;
   gpTMax = 1.99;
   h->SetPtRange(gpTMin, gpTMax);
   hMixed->SetPtRange(gpTMin, gpTMax);
  
   TH2* hist1 = 0;
   
-  c->cd(3);
-  GetDistAndFlow(h, 0, &hist1,  0, 6, 0, 20, 1.01, 4.99, 1);
-  hist1->DrawCopy("SURF1");
-
-  c->cd(6);
-  GetDistAndFlow(hMixed, 0, &hist1,  0, 6, 0, 20, 1.01, 4.99, 1);
-  hist1->DrawCopy("SURF1");
+  //h->GetUEHist(2)->SetSkipScaleMixedEvent(kTRUE);
   
-  c->cd(9);
-  GetSumOfRatios(h, hMixed, &hist1, 6, 0, 20, 1.01, 4.99, kTRUE); 
-  if (hist1)
-    hist1->DrawCopy("SURF1");  
+  if (h->GetUEHist(2)->GetTrackHist(0)->GetNVar() > 5)
+  {
+    Int_t step = 8;
+    
+    c->cd(3);
+    GetDistAndFlow(h, 0, &hist1,  0, step, 0, 20, 2.01, 2.99, 1);
+    hist1->DrawCopy("SURF1");
+
+    c->cd(6);
+    GetDistAndFlow(hMixed, 0, &hist1,  0, step, 0, 20, 2.01, 2.99, 1);
+    hist1->DrawCopy("SURF1");
+    
+    c->cd(9);
+    GetSumOfRatios(h, hMixed, &hist1, step, 0, 20, 2.01, 2.99, kTRUE); 
+    if (hist1)
+      hist1->DrawCopy("SURF1");  
+  }
 
   c->SaveAs(Form("qa/%s", c->GetTitle()));
 }
@@ -12647,14 +12833,23 @@ void CompareMixedEvent(const char* fileName)
 
 void FillParentTHnSparse(const char* fileName, Bool_t reduce = kFALSE, const char* tag = "")
 {
-  if (TString(fileName).BeginsWith("alien:"))
+  TString fileNameStr(fileName);
+  
+  if (fileNameStr.BeginsWith("/alice"))
+  {
+    fileNameStr = Form("alien://%s", fileNameStr.Data());
+    if (fileNameStr.EndsWith("merge"))
+      fileNameStr += "/AnalysisResults.root";
+  }
+  
+  if (fileNameStr.BeginsWith("alien:"))
     TGrid::Connect("alien:");
   
   loadlibs();
 
   TList* list = 0;
   
-  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName, &list, kFALSE, tag);
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileNameStr, &list, kFALSE, tag);
   Printf("We have %d axes", ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0)->GetNVar()));
   
   if (reduce)
@@ -12662,15 +12857,15 @@ void FillParentTHnSparse(const char* fileName, Bool_t reduce = kFALSE, const cha
   ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0))->FillParent();
   ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0))->DeleteContainers();
   
-  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE, tag);
+  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileNameStr, 0, kTRUE, tag);
   if (reduce)
     ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->ReduceAxis();
   ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->FillParent();
   ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->DeleteContainers();
   
-  TString newFileName(fileName);
+  TString newFileName(fileNameStr);
 
-  if (TString(fileName).BeginsWith("alien:"))
+  if (fileNameStr.BeginsWith("alien:"))
     newFileName = gSystem->BaseName(newFileName);
   
   newFileName.ReplaceAll(".root", "");
@@ -12849,8 +13044,8 @@ void PlotCorrections(const char* fileName, const char* tag = "")
   c = new TCanvas("c", "c", 1200, 800);
   c->Divide(3, 3);
 
-  c2 = new TCanvas("c2", "c2", 800, 600);
-  c2->Divide(2, 2);
+  c2 = new TCanvas("c2", "c2", 800, 900);
+  c2->Divide(2, 3);
 
   h->SetEtaRange(-0.89, 0.89);
 //   h->SetEtaRange(-1.19, 1.19);
@@ -12902,6 +13097,20 @@ void PlotCorrections(const char* fileName, const char* tag = "")
   
   h->SetZVtxRange(0, -1);
 
+  for (Int_t i=0; i<4; i++)
+  {
+    c2->cd(5);
+    h->SetPartSpecies(i);
+    proj = h->GetUEHist(2)->GetTrackingEfficiency(1);
+    proj->GetYaxis()->SetTitle("tracking efficiency");
+    proj->SetTitle(""); proj->SetStats(0);
+//     proj->GetXaxis()->SetTitle("#eta");
+    proj->SetLineColor(i+1);
+    proj->DrawClone((i == 0) ? "" : "SAME");
+  }
+  
+  return;
+  
 /*  c->cd(9);
   h->GetUEHist(2)->GetTrackingContamination()->Draw("COLZ");*/
   
@@ -12910,8 +13119,6 @@ void PlotCorrections(const char* fileName, const char* tag = "")
   new TCanvas;
   proj2->Draw();
   
-  return;
-  
   new TCanvas;
   hist = h->GetUEHist(2)->GetCorrelatedContamination();
 //   if (hist->GetEntries() > 0)
@@ -13098,7 +13305,7 @@ void SaveEfficiencyCorrection(const char* fileName, const char* tag = "", Bool_t
   
   Printf("%f", effCorr->GetEntries());
   
-  TObjString tag2(Form("corrections from file %s", fileName));
+  TObjString tag2(Form("corrections from file %s with tag %s", fileName, tag));
 
   file = TFile::Open("correction.root", "RECREATE");
   effCorr->Write();
@@ -13955,6 +14162,14 @@ void CondenseCentrality(const char* fileName, Float_t targetValue, Int_t step =
   h->GetUEHist(2)->CondenseBin(step, 3, 1, targetValue, from, to);
   hMixed->GetUEHist(2)->CondenseBin(step, 3, 1, targetValue, from, to);
   
+  Float_t events = h->GetCentralityDistribution()->Integral();
+  h->GetCentralityDistribution()->Reset();
+  h->GetCentralityDistribution()->Fill(targetValue, events);
+  
+  events = hMixed->GetCentralityDistribution()->Integral();
+  hMixed->GetCentralityDistribution()->Reset();
+  hMixed->GetCentralityDistribution()->Fill(targetValue, events);
+  
   TString newFileName(fileName);
   newFileName.ReplaceAll(".root", "");
   newFileName += "_condensed.root";
@@ -14680,3 +14895,16 @@ MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc, Int_t ntrd = 7)//ntrd is 7
   Float_t scale = (ntrd * 0.14638485 + (18 - ntrd) * 0.02406) / 18.;
   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),scale/0.03471));
 }
+
+void Browse(const char* fileName)
+{
+  loadlibs();
+
+  if (!gGrid && TString(fileName).BeginsWith("alien://"))
+    TGrid::Connect("alien://");
+  
+  TFile::Open(fileName);
+
+  new TBrowser;
+}
+
index ea99de7..511afec 100644 (file)
@@ -36,10 +36,10 @@ void DrawLatex(Float_t x, Float_t y, Int_t color, const char* text, Float_t text
   latex->Draw();
 }
 
-// 0    1    2    3     4     5     6      7      8          9          10         11         12         13         14          15
-// norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,moment3phi,moment3eta,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
-//            16   17   18   19    20    21    22     23     24         25         26         27         28         29         30          31         
-// second fit:norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,moment3phi,moment3eta,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
+// 0    1    2    3     4     5     6      7      8          9          10     11     12         13         14          15
+// norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,phirms,etarms,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
+//            16   17   18   19    20    21    22     23     24         25         26     27     28         29         30          31         
+// second fit:norm,dphi,deta,norm2,dphi2,deta2,chi2_1,chi2_2,moment2phi,moment2eta,phirms,etarms,moment4phi,moment4eta,kurtosisphi,kurtosiseta,
 // 32     33              34      
 // IAAFit,Yield(integral),IAAHist,
 //        35      36       37    38    39      40       41    42    43       44
@@ -479,7 +479,7 @@ void SubtractEtaGap1D(TH1* projPhi, TH1* projPhiSubtractPositive, TH1* projPhiSu
 void CalculateMomentsKurtosis(Float_t momentCalcLimit, TH1* proj, Int_t graphIDStart, Int_t histId, Float_t x, Float_t xE)
 {
   //   momentCalcLimit = 0.6;
-  for (Int_t n=2; n <= 4; n++)
+  for (Int_t n=2; n <= 4; n+=2) // skipping moment 3
   {
     Float_t moment = 0;
     Float_t sum = 0;
@@ -506,17 +506,101 @@ void CalculateMomentsKurtosis(Float_t momentCalcLimit, TH1* proj, Int_t graphIDS
   AddPoint(graphs[graphIDStart+6][histId], x, proj->GetKurtosis(1), xE, proj->GetKurtosis(11));
 }
 
+void AddRMS(TGraphErrors* graph, Float_t x, Float_t xE, TF1* func, TMatrixDSym& cov, Int_t sigma1, Int_t sigma2, Int_t weight)
+{
+  Float_t a = TMath::Abs(func->GetParameter(sigma1));
+  Float_t b = TMath::Abs(func->GetParameter(sigma2));
+  Float_t w = func->GetParameter(weight);
+
+  Float_t rms = a * w + b * (1.0 - w);
+
+  /*
+  func->Print(); Printf("%f %f %f", a, b, w); cov.Print();
+  Printf("%f %f %f %f %f %f", TMath::Power(w * func->GetParError(sigma1), 2),
+    TMath::Power((1.0 - w) * func->GetParError(sigma2), 2),
+    TMath::Power((a - b) * func->GetParError(weight), 2),
+    2 * (a - b) * w * cov(weight, sigma1),
+    2 * (a - b) * (1.0 - w) * cov(weight, sigma2),
+    2 * w * (1.0 - w) * cov(sigma1, sigma2));
+  */
+  
+  Float_t sigma_rms = 
+    TMath::Power(w * func->GetParError(sigma1), 2) + 
+    TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
+    TMath::Power((a - b) * func->GetParError(weight), 2) + 
+    2 * (a - b) * w * cov(weight, sigma1) +
+    2 * (a - b) * (1.0 - w) * cov(weight, sigma2) +
+    2 * w * (1.0 - w) * cov(sigma1, sigma2);
+    
+   if (sigma_rms < 0)
+   {
+     Printf("WARNING: error is negative!");
+     
+     sigma_rms = TMath::Power(w * func->GetParError(sigma1), 2) + 
+     TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
+     TMath::Power((a - b) * func->GetParError(weight), 2);
+   }
+    
+  sigma_rms = TMath::Sqrt(sigma_rms);
+  
+  Printf("RMS: %f +- %f", rms, sigma_rms);
+  AddPoint(graph, x, rms, xE, sigma_rms);
+}
+
+void AddKurtosis(TGraphErrors* graph, Float_t x, Float_t xE, TF1* func, TMatrixDSym& cov, Int_t sigma1, Int_t sigma2, Int_t weight)
+{
+  Float_t a = TMath::Abs(func->GetParameter(sigma1));
+  Float_t b = TMath::Abs(func->GetParameter(sigma2));
+  Float_t w = func->GetParameter(weight);
+
+  Float_t rms = a * w + b * (1.0 - w);
+  Float_t kurtosis = 3 * (w * TMath::Power(a, 4) + (1.0 - w) * TMath::Power(b, 4)) / rms / rms - 3;
+  
+  Float_t der_a = 12*a*a*a*w / TMath::Power(b*b+(a*a-b*b)*w, 2) - (12*a*w*(b*b*b*b+(a*a*a*a-b*b*b*b)*w))/TMath::Power(b*b+(a*a-b*b)*w, 3);
+  Float_t der_b = (3*(4*b*b*b-4*b*b*b*w))/TMath::Power((b*b+(a*a-b*b)*w), 2) - (6*(2*b-2*b*w)*(b*b*b*b+(a*a*a*a-b*b*b*b)* w))/TMath::Power(b*b+(a*a-b*b)*w, 3);
+  Float_t der_w = (3*(a*a*a*a-b*b*b*b))/TMath::Power(b*b+(a*a-b*b)*w, 2) - (6*(a*a-b*b)*(b*b*b*b+(a*a*a*a-b*b*b*b)*w))/TMath::Power(b*b+(a*a-b*b)*w,3);
+  
+  Float_t sigma_kurtosis = 
+    TMath::Power(der_a * func->GetParError(sigma1), 2) + 
+    TMath::Power(der_b * func->GetParError(sigma2), 2) + 
+    TMath::Power(der_w * func->GetParError(weight), 2) + 
+    2 * der_w * der_a * cov(weight, sigma1) +
+    2 * der_w * der_b * cov(weight, sigma2) +
+    2 * der_a * der_b * cov(sigma1, sigma2);
+    
+//   if (sigma_rms < 0)
+//   {
+//     Printf("WARNING: error is negative!");
+//     
+//     sigma_rms = TMath::Power(w * func->GetParError(sigma1), 2) + 
+//     TMath::Power((1.0 - w) * func->GetParError(sigma2), 2) + 
+//     TMath::Power((a - b) * func->GetParError(weight), 2);
+//   }
+    
+  sigma_kurtosis = TMath::Sqrt(sigma_kurtosis);
+  
+  Printf("Kurtosis: %f +- %f", kurtosis, sigma_kurtosis);
+  AddPoint(graph, x, kurtosis, xE, sigma_kurtosis);
+}
+
 Float_t kEtaLimit = 1.0;
 Float_t kOuterLimit = 1.59;
 
 
-Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int_t graphID, Float_t x, Float_t xE, Float_t yPosChi2, Bool_t quick, Int_t histId, Int_t limits, Bool_t twoTrack)
+Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int_t graphID, Float_t x, Float_t xE, Float_t yPosChi2, Bool_t quick, Int_t histId, Int_t limits, Bool_t /*twoTrack*/)
 {
   Bool_t success = kTRUE;
   Float_t etaLimit = kEtaLimit;
   Float_t outerLimit = kOuterLimit;
 
-  if (graphID >= 15)
+  // for pt,T < 3
+  if (graphID <= 10)
+    etaLimit = 1.4;
+  if (graphID <= 15)
+    etaLimit = 1.2;
+  
+  // for pT,a > 4
+  if (graphID == 18 || graphID == 23 || graphID == 24)
   {
     etaLimit = 0.5;
     outerLimit = 0.99;
@@ -529,6 +613,28 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
     etaFitUpperLimit = 0.6;
     initSigma = 0.4;
   }
+  if ((graphID == 10 && histId == 4) || (graphID == 5 && histId == 0))
+    etaFitUpperLimit += 0.1;
+  
+  Printf("Limits: etaLimit = %f outerLimit = %f sigmaFitLimit = %f etaFitUpperLimit = %f initSigma = %f", etaLimit, outerLimit, sigmaFitLimit, etaFitUpperLimit, initSigma);
+
+  // set errors large for bins around (0,0)
+  if (1 && graphID <= 12)
+  {
+    Float_t exclusionRegion = 0.19;
+    if (graphID == 0)
+      exclusionRegion = 0.29;
+    if (graphID == 10 || graphID == 11 || graphID == 12)
+      exclusionRegion = 0.075;
+    
+    Printf("NOTE : Skipping bins at (0, 0)");
+    for (Int_t binx = hist->GetXaxis()->FindBin(-exclusionRegion); binx <= hist->GetXaxis()->FindBin(exclusionRegion); binx++)
+      for (Int_t biny = hist->GetYaxis()->FindBin(-exclusionRegion); biny <= hist->GetYaxis()->FindBin(exclusionRegion); biny++)
+      {
+//     hist->SetBinContent(binx, biny,  0);
+       hist->SetBinError(binx, biny,  1e5);
+      }
+  }
 
   hist->GetYaxis()->SetRangeUser(-outerLimit+0.01, outerLimit-0.01);
   hist->GetXaxis()->SetRangeUser(-TMath::Pi() / 2 + 0.01, TMath::Pi() * 0.5 - 0.01);
@@ -574,10 +680,10 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
     //   func->SetParameters(1, 0.3, 0.3, 0.25, initSigma, initSigma);
     func->SetParLimits(0, 0, 10);
     func->SetParLimits(1, sigmaFitLimit, 0.6);
-    func->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit + (((graphID == 5 && histId == 4) || (graphID == 0 && histId == 0)) ? 0.1 : 0));
-    func->SetParLimits(3, 0.1, 0.9);
-    func->SetParLimits(4, sigmaFitLimit, 0.6);
-    func->SetParLimits(5, sigmaFitLimit, etaFitUpperLimit + (((graphID == 5 && histId == 4) || (graphID == 0 && histId == 0)) ? 0.1 : 0));
+    func->SetParLimits(2, sigmaFitLimit, etaFitUpperLimit);
+    func->SetParLimits(3, 0.05, 0.95);
+    func->SetParLimits(4, sigmaFitLimit, 0.601);
+    func->SetParLimits(5, sigmaFitLimit, etaFitUpperLimit);
     for (Int_t i=6; i<bins+6; i++)
       func->FixParameter(i, func->GetParameter(i));
     hist->GetYaxis()->SetRangeUser(-etaLimit+0.01, etaLimit-0.01);
@@ -602,9 +708,10 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
     if (fitResult != 0)
       success = kFALSE;
   }
-  Int_t fitResult = hist->Fit(func, "0", "");
-  Printf("Fit result: %d; Chi2/ndf: %f/%d", fitResult, func->GetChisquare(), func->GetNDF());
-  if (fitResult != 0)
+  TFitResultPtr fitResultPtr = hist->Fit(func, "S0", "");
+  TMatrixDSym cov = fitResultPtr->GetCovarianceMatrix();
+  Printf("Fit result: %d; Chi2/ndf: %f/%d", (Int_t) fitResultPtr, func->GetChisquare(), func->GetNDF());
+  if (fitResultPtr != 0)
     success = kFALSE;
 
   Printf("Trying 1 Gaussian...");
@@ -620,8 +727,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   func_clone->FixParameter(3, 1);
   func_clone->FixParameter(4, sigmaFitLimit);
   func_clone->FixParameter(5, sigmaFitLimit);
-  fitResult = hist->Fit(func_clone, "0R", "");
-  Printf("Fit result: %d", fitResult);
+  TFitResultPtr fitResultPtr2 = hist->Fit(func_clone, "S0R", "");
+  Printf("Fit result: %d", (Int_t) fitResultPtr2);
   
   // if both parameters are within 1%, refit with 1 Gaussian only
   if (TMath::Abs(1.0 - func->GetParameter(1) / func->GetParameter(4)) < 0.01 && TMath::Abs(1.0 - func->GetParameter(2) / func->GetParameter(5)) < 0.01)
@@ -629,8 +736,9 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
     Printf("Parameters within 1%%. Using result from 1 Gaussian...");
     
     func = func_clone;
+    cov = fitResultPtr2->GetCovarianceMatrix();
     
-    if (fitResult != 0)
+    if (fitResultPtr2 != 0)
       success = kFALSE;
   }
   else if (func_clone->GetChisquare() * 0.99 < func->GetChisquare())
@@ -638,8 +746,9 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
     Printf("1 Gaussian fit has a at maximum 1%% (%f %f) larger chi2. Using results from 1 Gaussian...", func_clone->GetChisquare(), func->GetChisquare());
     
     func = func_clone;
+    cov = fitResultPtr2->GetCovarianceMatrix();
     
-    if (fitResult != 0)
+    if (fitResultPtr2 != 0)
       success = kFALSE;
   } 
   
@@ -665,6 +774,12 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   else
     AddPoint(graphs[3+16][graphID], x, 1.0 - func->GetParameter(3), xE, func->GetParError(3));    
   
+  //rms
+  AddRMS(graphs[10+16][graphID], x, xE, func, cov, 1, 4, 3);
+  AddRMS(graphs[11+16][graphID], x, xE, func, cov, 1+1, 4+1, 3);
+  AddKurtosis(graphs[14+16][graphID], x, xE, func, cov, 1, 4, 3);
+  AddKurtosis(graphs[15+16][graphID], x, xE, func, cov, 1+1, 4+1, 3);
+  
   //   hist->Fit(func, "0RI", "");
 
   //   width[0]->SetPoint(width[0]->GetN(), x, TMath::Abs(func->GetParameter(1)));
@@ -696,7 +811,7 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   residuals->SetMaximum((max - min) / 2);
   residuals->Draw("SURF1");
 
-  Float_t chi2 = 0;
+  Double_t chi2 = 0;
   Int_t ndf = 0;
   for (Int_t i=hist->GetXaxis()->FindBin(-0.8); i<=hist->GetXaxis()->FindBin(0.8); i++)
     for (Int_t j=hist->GetYaxis()->FindBin(-etaLimit); j<=hist->GetYaxis()->FindBin(etaLimit); j++)
@@ -733,8 +848,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   funcHist = (TH2*) hist->Clone("funcHistb");
   funcHist->Reset();
   funcHist->Add(funcClone);
-  funcHist->SetMinimum(-(max - min) / 2);
-  funcHist->SetMaximum((max - min) / 2);
+  funcHist->SetMinimum(0);
+  funcHist->SetMaximum(max - min);
   funcHist->Draw("SURF1");
   
   // eta gap subtraction
@@ -742,8 +857,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   func->SetParameter(0, 0);
   TH2* subtractFlow = (TH2*) hist->Clone("subtractFlow");
   subtractFlow->Add(func, -1);
-  subtractFlow->SetMinimum(-(max - min) / 2);
-  subtractFlow->SetMaximum((max - min) / 2);
+  subtractFlow->SetMinimum(0);
+  subtractFlow->SetMaximum(max - min);
   subtractFlow->DrawCopy("SURF1");
   sumSummary->SetPoint(sumSummary->GetN(), sumSummary->GetN(), subtractFlow->Integral());
 /*
@@ -762,7 +877,7 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
 */
 
   Float_t momentFitLimit = 0.8 - 1e-4;
-  if (graphID >= 15)
+  if (graphID >= 20)
     momentFitLimit = 0.4 - 1e-4;
 
   TH1* projx3 = hist->ProjectionX(Form("%s_projx3", hist->GetName()), hist->GetYaxis()->FindBin(-etaLimit+0.01), hist->GetYaxis()->FindBin(etaLimit-0.01));
@@ -848,8 +963,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   nBins = hist->GetXaxis()->FindBin(momentFitLimit-0.01) - hist->GetXaxis()->FindBin(-momentFitLimit+0.01)+1;
   projy2->Scale(1.0/nBins);
   
-  CalculateMomentsKurtosis(momentFitLimit, projx2, 8, graphID, x, xE);
-  CalculateMomentsKurtosis(momentFitLimit, projy2, 9, graphID, x, xE);
+//   CalculateMomentsKurtosis(momentFitLimit, projx2, 8, graphID, x, xE);
+//   CalculateMomentsKurtosis(momentFitLimit, projy2, 9, graphID, x, xE);
 
 //     return success;
   
@@ -859,8 +974,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   TH1* projy1 = subtractFlow->ProjectionY(Form("%s_projy1", hist->GetName()), hist->GetXaxis()->FindBin(-0.79), hist->GetXaxis()->FindBin(0.79));
   projy1->GetXaxis()->SetRangeUser(-momentFitLimit, momentFitLimit);
 
-  CalculateMomentsKurtosis(momentFitLimit, projx1, 8+16, graphID, x, xE);
-  CalculateMomentsKurtosis(momentFitLimit, projy1, 9+16, graphID, x, xE);
+//   CalculateMomentsKurtosis(momentFitLimit, projx1, 8+16, graphID, x, xE);
+//   CalculateMomentsKurtosis(momentFitLimit, projy1, 9+16, graphID, x, xE);
 
 //   TF1* twoGauss = new TF1("twoGauss", "gaus(0)+gaus(3)", -2, 2);
 //   twoGauss->SetParameters(1, 0, 0.3, 1, 0, 0.6);
@@ -903,19 +1018,6 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   // 1d fit (lots of params)
   AddPoint(etaWidthSummary, etaWidthSummary->GetN(), projy1->GetFunction("gaus")->GetParameter(2), 0, projy1->GetFunction("gaus")->GetParError(2));
 
-  // set errors large for bins potentially affected by two-track effects
-  if (0 && twoTrack)
-  {
-    Printf("NOTE : Skipping bins at (0, 0)");
-//     for (Int_t binx = hist->GetXaxis()->FindBin(-0.25); binx <= hist->GetXaxis()->FindBin(0.25); binx++)
-    for (Int_t binx = hist->GetXaxis()->FindBin(-0.01); binx <= hist->GetXaxis()->FindBin(0.01); binx++)
-      for (Int_t biny = hist->GetYaxis()->FindBin(-0.01); biny <= hist->GetYaxis()->FindBin(0.01); biny++)
-      {
-//     hist->SetBinContent(binx, biny,  0);
-       hist->SetBinError(binx, biny,  1e5);
-      }
-  }
-
   Float_t etaFitLimit = outerLimit;
 
   Bool_t oneGaussian = kTRUE;
@@ -955,8 +1057,8 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   }
   canvas->cd(canvasPos++);
   func4phi->SetLineColor(4);
-  TFitResultPtr fitResultPtr = projx3->Fit(func4phi, "SIR", "SAME");
-  TMatrixDSym cov = fitResultPtr->GetCovarianceMatrix();
+  fitResultPtr = projx3->Fit(func4phi, "SIR", "SAME");
+  TMatrixDSym cov2 = fitResultPtr->GetCovarianceMatrix();
 
   projx3->Draw("SAME");
   projx3->SetLineColor(4);
@@ -995,7 +1097,7 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   {
     Float_t vector[3];
     for (Int_t i=2; i<5;i++)
-      vector[i-2] = cov[i][2]*func4phi->GetParameter(4) + cov[i][3]*(1-func4phi->GetParameter(4))+cov[i][4]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3)));
+      vector[i-2] = cov2[i][2]*func4phi->GetParameter(4) + cov2[i][3]*(1-func4phi->GetParameter(4))+cov2[i][4]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3)));
 
     Float_t sigma = TMath::Sqrt(vector[0]*func4phi->GetParameter(4) + vector[1]*(1-func4phi->GetParameter(4))+ vector[2]*(TMath::Abs(func4phi->GetParameter(2))-TMath::Abs(func4phi->GetParameter(3))));
     Float_t rms = TMath::Abs(func4phi->GetParameter(2))*func4phi->GetParameter(4)+TMath::Abs(func4phi->GetParameter(3))*(1-func4phi->GetParameter(4));
@@ -1029,7 +1131,7 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   canvas->cd(canvasPos++);
   func4eta->SetLineColor(4);
   fitResultPtr = projy3->Fit(func4eta, "SIR", "SAME");
-  cov = fitResultPtr->GetCovarianceMatrix();
+  cov2 = fitResultPtr->GetCovarianceMatrix();
 
   projy3->Draw("SAME");
   projy3->SetLineColor(4);
@@ -1044,6 +1146,9 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   for (Int_t i=projy3->FindBin(-momentFitLimit); i<=projy3->FindBin(momentFitLimit);i++)
   {
 //    Float_t temp = chi2;
+    if (projy3->GetBinError(i) <= 0)
+      continue;
+    
     chi2 += TMath::Power((projy3->GetBinContent(i)-func4eta->Eval(projy3->GetBinCenter(i)))/projy3->GetBinError(i),2);
 //    cerr << i << "\t" << projy3->GetBinCenter(i) << "\t" << chi2-temp << endl;
     ndf++;
@@ -1065,7 +1170,7 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   {
     Float_t vector[3];
     for (Int_t i=2; i<5;i++)
-      vector[i-2] = cov[i][2]*func4eta->GetParameter(4) + cov[i][3]*(1-func4eta->GetParameter(4))+cov[i][4]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3)));
+      vector[i-2] = cov2[i][2]*func4eta->GetParameter(4) + cov2[i][3]*(1-func4eta->GetParameter(4))+cov2[i][4]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3)));
 
     Float_t sigma = TMath::Sqrt(vector[0]*func4eta->GetParameter(4) + vector[1]*(1-func4eta->GetParameter(4))+ vector[2]*(TMath::Abs(func4eta->GetParameter(2))-TMath::Abs(func4eta->GetParameter(3))));
     Float_t rms = TMath::Abs(func4eta->GetParameter(2))*func4eta->GetParameter(4)+TMath::Abs(func4eta->GetParameter(3))*(1-func4eta->GetParameter(4));
@@ -1104,14 +1209,17 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
 //   Float_t etaFitLimit = 0.5;
   TF2* func3 = new TF2("func3", "[0]+[1]*([4]/TMath::TwoPi()/[2]/[3]*exp(-0.5*((x/[2])**2+(y/[3])**2))+(1-[4])/TMath::TwoPi()/[5]/[6]*exp(-0.5*((x/[5])**2+(y/[6])**2)))", -0.5 * TMath::Pi(), 0.5 * TMath::Pi(), -etaFitLimit, etaFitLimit);
 
-  func3->SetParLimits(4, 0.1, 0.9);
+  func3->SetParameters(0, hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)), 0.3, 0.3, 0.25, initSigma, initSigma);
+
+  func3->FixParameter(0, 0);
+
+  func3->SetParLimits(4, 0.05, 0.95);
   func3->SetParLimits(1, 0, 10);
-  func3->SetParLimits(2, sigmaFitLimit, 0.7);
+  func3->SetParLimits(2, sigmaFitLimit, 0.6);
   func3->SetParLimits(3, sigmaFitLimit, etaFitUpperLimit);
-  func3->SetParLimits(5, sigmaFitLimit, 0.7);
+  func3->SetParLimits(5, sigmaFitLimit, 0.601);
   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
-  
-  func3->FixParameter(0, 0);
+
 /*
   sigmaFitLimit = 0.09;
   etaFitUpperLimit = 0.12;
@@ -1123,29 +1231,16 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   func3->SetParLimits(6, sigmaFitLimit, etaFitUpperLimit);
   func3->FixParameter(0, 0);
 */
-  func3->SetParameters(0, hist->GetBinContent(hist->GetXaxis()->FindBin(0.0), hist->GetYaxis()->FindBin(0.0)) - mean, 0.3, 0.3, 0.25, initSigma, initSigma);
 //   for (Int_t i=0; i<6; i++)
 //     func3->SetParameter(i+1, func->GetParameter(i));
 
-  if (0 && histId == 0)
-  {
-    // central --> go to 1 Gaussian only
-    func3->SetParLimits(4, 1, 1);
-    func3->FixParameter(4, 1);
-    func3->FixParameter(5, sigmaFitLimit);
-    func3->FixParameter(6, sigmaFitLimit);
-  }
-  
-  // set errors 20% of the value for bins potentially affected by two-track effects
-//   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(0.0001)));
-//   hist->SetBinError(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(0.0001),  hist->GetYaxis()->FindBin(-0.0001)));
-//   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001),  0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(0.0001)));
-//   hist->SetBinError(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001), 0.1 * hist->GetBinContent(hist->GetXaxis()->FindBin(-0.0001), hist->GetYaxis()->FindBin(-0.0001)));
-
-  fitResult = hist->Fit(func3, "0R", "");
-  Printf("Fit result: %d", fitResult);
-  if (fitResult != 0)
+  fitResultPtr = hist->Fit(func3, "0RS", "");
+  Printf("Fit result: %d", (Int_t) fitResultPtr);
+//   return 0;
+  if (fitResultPtr != 0)
     success = kFALSE;
+  TMatrixDSym cov3 = fitResultPtr->GetCovarianceMatrix();
+  cov3.Print();
   
   Printf("Testing 1 Gaussian...");
 
@@ -1156,25 +1251,28 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   func3_clone->FixParameter(5, sigmaFitLimit);
   func3_clone->FixParameter(6, sigmaFitLimit);
   
-  fitResult = hist->Fit(func3_clone, "0R", "");
-  Printf("Fit result: %d", fitResult);
+  fitResultPtr2 = hist->Fit(func3_clone, "0RS", "");
+  Printf("Fit result: %d", (Int_t) fitResultPtr2);
   
   // if both parameters were within 1%, use 1 Gaussian parameters
   if (TMath::Abs(1.0 - func3->GetParameter(2) / func3->GetParameter(5)) < 0.01 && TMath::Abs(1.0 - func3->GetParameter(3) / func3->GetParameter(6)) < 0.01)
   {
     Printf("Parameters within 1%%. Using results from 1 Gaussian...");
     
-    if (fitResult != 0)
+    if (fitResultPtr2 != 0)
       success = kFALSE;
     func3 = func3_clone;
+    cov3 = fitResultPtr2->GetCovarianceMatrix();
+//     cov3.Print();
   }
   else if (func3_clone->GetChisquare() * 0.99 < func3->GetChisquare())
   {
-    Printf("1 Gaussian fit has a at maximum 1%% (%f %f) larger chi2. Using results from 1 Gaussian...", func3_clone->GetChisquare(), func3->GetChisquare());
+    Printf("1 Gaussian fit has a at maximum 1%% (%f/%d %f/%d) larger chi2. Using results from 1 Gaussian...", func3_clone->GetChisquare(), func3_clone->GetNDF(), func3->GetChisquare(), func3_clone->GetNDF());
     
-    if (fitResult != 0)
+    if (fitResultPtr2 != 0)
       success = kFALSE;
     func3 = func3_clone;
+    cov3 = fitResultPtr2->GetCovarianceMatrix();
   }
   
 //   hist->Fit(func3, "I0R", "");
@@ -1207,6 +1305,9 @@ Bool_t FitDeltaPhi2DOneFunction(TH2* hist, TCanvas* canvas, Int_t canvasPos, Int
   else
     AddPoint(graphs[3][graphID], x, 1.0 - func3->GetParameter(4), xE, func3->GetParError(4));
 
+  AddRMS(graphs[10][graphID], x, xE, func3, cov3, 2, 5, 4);
+  AddRMS(graphs[11][graphID], x, xE, func3, cov3, 2+1, 5+1, 4);
+  
   TH2* funcHist3 = (TH2*) hist->Clone("funcHist3");
   funcHist3->Reset();
   funcHist3->Add(func3);
@@ -1305,11 +1406,15 @@ void AnalyzeDeltaPhiEtaGap2D(const char* fileName, const char* outputFileName =
        continue;
       }
       
-      for (Int_t histId = 0; histId < NHists; histId++)
+      for (Int_t histId = 0; histId < 8; histId++)
       {
 //     if (i != 1 || j != 2)
 //       continue;
-    
+
+       // skip 4.0 < p_{T,trig} < 8.0 - 6.00 < p_{T,assoc} < 8.00 above 40% and pp
+       if (i == 2 && j == 5 && (histId == 1 || histId == 2 || histId == 5))
+         continue;
+       
        hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
        if (!hist1)
         {
@@ -1344,8 +1449,8 @@ void AnalyzeDeltaPhiEtaGap2D(const char* fileName, const char* outputFileName =
          for (Int_t k=0; k<NGraphs; k++)
            graphs[k][graphID]->SetTitle(hist1->GetTitle());
        
-       Float_t centralityAxisMapping[] = { 5, 65, 100, 30, 15, 50 };
-       Float_t centralityAxisMappingE[] = { 5, 5, 0, 10, 5, 10 };
+       Float_t centralityAxisMapping[] = { 5, 70, 100, 30, 15, 50 };
+       Float_t centralityAxisMappingE[] = { 5, 10, 0, 10, 5, 10 };
 
        Bool_t success = FitDeltaPhi2DOneFunction((TH2*) hist1, canvas, 1, graphID, centralityAxisMapping[histId], centralityAxisMappingE[histId], 0.9, kTRUE, histId, (i > 0) ? 1 : 0, (j <= 2 && histId != 2));
        if (!success)
@@ -1434,7 +1539,7 @@ void AnalyzeDeltaPhiEtaGap2DExample(const char* fileName, Int_t i, Int_t j, Int_
   paveText->Draw();
   
   c->cd(1);
-  DrawALICELogo(kTRUE, 0.2, 0.7, 0.4, 0.9);
+//   DrawALICELogo(kTRUE, 0.2, 0.7, 0.4, 0.9);
   
 //   return;
     
@@ -1592,16 +1697,20 @@ void CalculateRMS(Int_t nHists, TGraphErrors** graph, TGraphErrors** graph2, TGr
       graph[histId]->GetY()[i] = rms;
       // error**2 = weight**2 * error_sigma**2 + (weight-1)**2 * error_sigma2**2 + (sigma1 + sigma2)**2 * error_weight**2
       // TODO this neglects some correlations
-      graph[histId]->GetEY()[i] = TMath::Sqrt(TMath::Power(weight[histId]->GetY()[i] * graph[histId]->GetEY()[i], 2) + 
-       TMath::Power((weight[histId]->GetY()[i] - 1) * graph2[histId]->GetEY()[i], 2) + TMath::Power((graph[histId]->GetY()[i] + graph2[histId]->GetY()[i]) * weight[histId]->GetEY()[i], 2));
+      graph[histId]->GetEY()[i] = TMath::Sqrt(
+       TMath::Power(weight[histId]->GetY()[i] * graph[histId]->GetEY()[i], 2) + 
+       TMath::Power((weight[histId]->GetY()[i] - 1) * graph2[histId]->GetEY()[i], 2) + 
+       TMath::Power((graph[histId]->GetY()[i] - graph2[histId]->GetY()[i]) * weight[histId]->GetEY()[i], 2)
+      );
     }
   }
 }
 
 Bool_t SkipGraph(Int_t i)
 {
+  return (i == 14);
 //  return (i == 1 || i == 7 || i == 9 || i == 13 || i == 14);
-  return (i == 7 || i == 9 || i == 13 || i == 14); //For the STAR comparison
+//   return (i == 7 || i == 9 || i == 13 || i == 14); //For the STAR comparison
 }
 
 Int_t skipGraphList[] =  { -1, -1, -1 };
@@ -1799,12 +1908,12 @@ const char* MCLabel = 0;
 
 void DrawCentrality(const char* canvasName, Int_t nHists, TGraphErrors** graph, Float_t min = 0, Float_t max = 0, const char* yLabel = "", TGraphErrors** systematicA = 0, TGraphErrors** systematicB = 0, TGraphErrors** graph2 = 0, TGraphErrors** systematic2 = 0, TGraphErrors** graph3 = 0, Int_t uncertaintyID = -1)
 {
-  Bool_t found = kTRUE;
+//   Bool_t found = kTRUE;
   TCanvas* c1 = 0;//(TCanvas*) gROOT->GetListOfCanvases()->FindObject(canvasName);
   if (!c1)
   {
     c1 = new TCanvas(canvasName, canvasName, 800, 600);
-    found = kFALSE;
+//     found = kFALSE;
   }
   c1->cd();
   
@@ -1967,8 +2076,8 @@ void DrawCentrality(const char* canvasName, Int_t nHists, TGraphErrors** graph,
     DrawALICELogo(kTRUE, 0.41, 0.75, 0.54, 0.95);
   
   gSystem->mkdir(fgFolder, kTRUE);
-  c1->SaveAs(Form("%s/%s.gif", fgFolder.Data(), canvasName));
-  c1->SaveAs(Form("%s/%s.eps", fgFolder.Data(), canvasName));
+  c1->SaveAs(Form("%s/%s.png", fgFolder.Data(), canvasName));
+//   c1->SaveAs(Form("%s/%s.eps", fgFolder.Data(), canvasName));
 }
 
 void CalculateRMSSigma(TGraphErrors*** graphsTmp = 0)
@@ -2046,13 +2155,17 @@ void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fil
 //    DrawCentrality("norm2phi(1D)", nHists, graphs[36], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
 //    DrawCentrality("norm2eta(1D)", nHists, graphs[40], 0, 1.5, "N (a.u.)", (graphsWingRemoved) ? graphsWingRemoved[0+offset] : 0);
 
+    
+    Bool_t uncertainties = kFALSE;
+
     CalculateRMSSigma();
     if (graphsWingRemoved)
       CalculateRMSSigma(graphsWingRemoved);
+//     DrawCentrality("phi_rms_centrality_old", nHists, graphs[1+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
+//     DrawCentrality("eta_rms_centrality_old", nHists, graphs[2+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
 
-    DrawCentrality("phi_rms_centrality", nHists, graphs[1+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, 1);
-
-    DrawCentrality("eta_rms_centrality", nHists, graphs[2+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0, 0, 0, 0, 0, 1);
+    DrawCentrality("phi_rms_centrality", nHists, graphs[10+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[10+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
+    DrawCentrality("eta_rms_centrality", nHists, graphs[11+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[11+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
 
 //     return;
 //   DrawCentrality("phi_rms_centrality(1D)", nHists, graphs[37], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[1+offset] : 0, 0, 0, 0, 0, 1);
@@ -2060,13 +2173,13 @@ void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fil
 //   DrawCentrality("eta_rms_centrality(1D)", nHists, graphs[41], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), (graphsWingRemoved) ? graphsWingRemoved[2+offset] : 0, 0, 0, 0, 0, 1);
 //   DrawCentrality("chi2_phi(1D)", nHists, graphs[43], 0, 150, "#chi^{2}/ndf");
 //   DrawCentrality("chi2_eta(1D)", nHists, graphs[44], 0, 150, "#chi^{2}/ndf");
-/*    
-    DrawCentrality("chi2_1", nHists, graphs[6+offset], 0.5, 5, "#chi^{2}/ndf (full region)");
-    DrawCentrality("chi2_2", nHists, graphs[7+offset], 0.5, 5, "#chi^{2}/ndf (peak region)");
+    DrawCentrality("chi2_1", nHists, graphs[6+offset], 0.5, 5, "#chi^{2}/ndf (full region)", (graphsWingRemoved) ? graphsWingRemoved[6+offset] : 0);
+    DrawCentrality("chi2_2", nHists, graphs[7+offset], 0.5, 5, "#chi^{2}/ndf (peak region)", (graphsWingRemoved) ? graphsWingRemoved[7+offset] : 0);
     
-    DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 1, "#sigma_{#Delta#varphi} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[8+offset] : 0, 0, 0, 0, 0, 1);
-    DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 1, "#sigma_{#Delta#eta}", (graphsWingRemoved) ? graphsWingRemoved[9+offset] : 0, 0, 0, 0, 0, 1);
+    DrawCentrality("sigma_phi", nHists, graphs[8+offset], 0, 1, "#sigma_{#Delta#varphi} (rad.)", (graphsWingRemoved) ? graphsWingRemoved[8+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
+    DrawCentrality("sigma_eta", nHists, graphs[9+offset], 0, 1, "#sigma_{#Delta#eta}", (graphsWingRemoved) ? graphsWingRemoved[9+offset] : 0, 0, 0, 0, 0, (uncertainties) ? 1 : -1);
 
+/*    
 //     DrawCentrality("moment3_phi", nHists, graphs[10+offset], -0.5, 0.5, "moment3 #varphi (rad.)", (graphsWingRemoved) ? graphsWingRemoved[10+offset] : 0, 0, 0, 0, 0, 0);
 //     DrawCentrality("moment3_eta", nHists, graphs[11+offset], -0.5, 0.5, "moment3 #eta", (graphsWingRemoved) ? graphsWingRemoved[11+offset] : 0, 0, 0, 0, 0, 0);
 */
@@ -2077,8 +2190,8 @@ void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fil
 */
   }
   
-//  DrawCentrality("kurtosisphi_centrality", 12, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", (graphsWingRemoved) ? graphsWingRemoved[14+offset] : 0, 0, 0, 0, 0, 2);
-//  DrawCentrality("kurtosiseta_centrality", 12, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", (graphsWingRemoved) ? graphsWingRemoved[15+offset] : 0, 0, 0, 0, 0, 2);
+ DrawCentrality("kurtosisphi_centrality", 12, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", (graphsWingRemoved) ? graphsWingRemoved[14+offset] : 0, 0, 0, 0, 0, 2);
+ DrawCentrality("kurtosiseta_centrality", 12, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", (graphsWingRemoved) ? graphsWingRemoved[15+offset] : 0, 0, 0, 0, 0, 2);
 }
 
 /*
@@ -2129,7 +2242,7 @@ void DrawResultsCentrality(const char* fileName = "graphs.root", const char* fil
 
 void MCComparison(const char* fileNameData, const char* fileNameWingRemoved, const char* fileNameHijing, const char* fileNameAMPT, Int_t offset = 0)
 {
-  Int_t nHists = 12; //NHists;
+  Int_t nHists = 14; //NHists;
 
   ReadGraphs(fileNameWingRemoved);
   CalculateRMSSigma();
@@ -2153,7 +2266,7 @@ void MCComparison(const char* fileNameData, const char* fileNameWingRemoved, con
 
   if (0)
   {
-    Int_t graphList[] = { 1, 2, 8, 9, 14, 15 };
+    Int_t graphList[] = { 10, 11, 8, 9, 14, 15 };
     for (Int_t i=0; i<6; i++)
     {
       graphs[graphList[i]+offset][5] = new TGraphErrors;
@@ -2163,17 +2276,17 @@ void MCComparison(const char* fileNameData, const char* fileNameWingRemoved, con
     }
   }
   
-  DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1+offset], 0.2, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), graphsWingRemoved[1+offset], 0, graphs1[1+offset], 0, graphs2[1+offset], 1);
+  DrawCentrality("phi_rms_centrality_mc", nHists, graphs[10+offset], 0, 0.8, Form("#sigma_{#Delta#varphi} (%s) (rad.)", fitLabel), graphsWingRemoved[1+offset], 0, graphs1[10+offset], 0, graphs2[10+offset], 1);
   
 //   return;
   
-  DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2+offset], 0.2, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), graphsWingRemoved[2+offset], 0, graphs1[2+offset], 0, graphs2[2+offset], 1);
+  DrawCentrality("eta_rms_centrality_mc", nHists, graphs[11+offset], 0, 0.8 + 0.4 / 16 * offset, Form("#sigma_{#Delta#eta} (%s)", fitLabel), graphsWingRemoved[11+offset], 0, graphs1[11+offset], 0, graphs2[2+offset], 1);
   
-  DrawCentrality("sigma_phi_mc", nHists, graphs[8+offset], 0.2, 0.6, "#sigma_{#Delta#varphi} (rad.)", graphsWingRemoved[8+offset], 0, graphs1[8+offset], 0, graphs2[8+offset], 1);
-  DrawCentrality("sigma_eta_mc", nHists, graphs[9+offset], 0.2, 0.6, "#sigma_{#Delta#eta}", graphsWingRemoved[9+offset], 0, graphs1[9+offset], 0, graphs2[9+offset], 1);
+  DrawCentrality("sigma_phi_mc", nHists, graphs[8+offset], 0, 0.6, "#sigma_{#Delta#varphi} (rad.)", graphsWingRemoved[8+offset], 0, graphs1[8+offset], 0, graphs2[8+offset], 1);
+  DrawCentrality("sigma_eta_mc", nHists, graphs[9+offset], 0, 0.6, "#sigma_{#Delta#eta}", graphsWingRemoved[9+offset], 0, graphs1[9+offset], 0, graphs2[9+offset], 1);
 
-  DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -1.5, 2.5, "Kurtosis #Delta#varphi", graphsWingRemoved[14+offset], 0, graphs1[14+offset], 0, graphs2[14+offset], 2);
-  DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -1.5, 2.5, "Kurtosis #Delta#eta", graphsWingRemoved[15+offset], 0, graphs1[15+offset], 0, graphs2[15+offset], 2);
+  DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14+offset], -1.5, 3.5, "Kurtosis #Delta#varphi", graphsWingRemoved[14+offset], 0, graphs1[14+offset], 0, graphs2[14+offset], 2);
+  DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15+offset], -1.5, 3.5, "Kurtosis #Delta#eta", graphsWingRemoved[15+offset], 0, graphs1[15+offset], 0, graphs2[15+offset], 2);
 
 /*  DrawCentrality("phi_rms_centrality_mc", nHists, graphs[1], 0, 0.8, "#sigma_{#varphi} (fit) (rad.)", graphs[1+16], graphsWingRemoved[1],  graphs1[1], 0, graphs2[1], 1);
   DrawCentrality("eta_rms_centrality_mc", nHists, graphs[2], 0, 0.8, "#sigma_{#eta} (fit)", graphs[2+16], graphsWingRemoved[2], graphs1[2], 0, graphs2[2], 1);
@@ -2208,7 +2321,7 @@ void ShowWingEffect(const char* fileNameData, const char* fileNameWingRemoved, I
 
 void ShowFitEffect(const char* fileNameData)
 {
-  Int_t nHists = 12; //NHists;
+  Int_t nHists = NHists;
 
   ReadGraphs(fileNameData);
   CalculateRMSSigma();
@@ -2221,6 +2334,9 @@ void ShowFitEffect(const char* fileNameData)
 
   DrawCentrality("kurtosisphi_centrality_mc", nHists, graphs[14], -2, 4, "Kurtosis #varphi", graphs[14+16]);
   DrawCentrality("kurtosiseta_centrality_mc", nHists, graphs[15], -2, 4, "Kurtosis #eta", graphs[15+16]);
+
+  DrawCentrality("chi2_1", nHists, graphs[6], 0.5, 5, "chi2_1", graphs[6+16]);
+  DrawCentrality("chi2_2", nHists, graphs[7], 0.5, 5, "chi2_2", graphs[7+16]);
 }
 
 void CompareSigmas(const char* fileNameData)
@@ -3049,9 +3165,70 @@ TH1* GetSystUnc(TH1* hist, Int_t i, Int_t j, Int_t histId, Int_t etaPhi)
   return systUnc;
 }
 
+void Draw2DExamples(const char* histFileName)
+{
+  Int_t is[] = { 0, 0, 1, 1, 2, 2, 3 };
+  Int_t js[] = { 1, 2, 2, 3, 4, 5, 6 };
+  Int_t cs[] = { 0, 3, 5, 1 };
+  Float_t outerLimit = kOuterLimit;
+
+  TFile::Open(histFileName);
+
+  TCanvas* canvas = new TCanvas("2d", "2d", 1800, 600);
+  canvas->Divide(7, 4);
+  Int_t cID = 1;
+
+  for (Int_t c=0; c<4; c++)
+  {
+    for (Int_t i=0; i<7; i++)
+    {
+      canvas->cd(cID++);
+
+      TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", is[i], js[i], cs[c]));
+      if (!hist)
+       continue;
+
+      // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+      hist->Scale(1.0 / hist->GetYaxis()->GetBinWidth(1));
+      hist->GetZaxis()->SetTitle(kCorrFuncTitle);
+      hist->GetZaxis()->SetTitleOffset(1.9);
+      
+      hist->Rebin2D(2, 2);
+      hist->Scale(0.25);
+  
+      TString label(hist->GetTitle());
+      label.ReplaceAll(".00", " GeV/c");
+      label.ReplaceAll(".0", " GeV/c");
+      TObjArray* objArray = label.Tokenize("-");
+      TPaveText* paveText = new TPaveText(0.52, 0.72, 0.95, 0.95, "BRNDC");
+      paveText->SetTextSize(0.035);
+      paveText->SetFillColor(0);
+      paveText->SetShadowColor(0);
+      paveText->AddText(objArray->At(0)->GetName());
+      paveText->AddText(objArray->At(1)->GetName());
+      if (objArray->GetEntries() == 4)
+       paveText->AddText(Form("Pb-Pb 2.76 TeV %s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
+      else
+       paveText->AddText(Form("%s 2.76 TeV", objArray->At(2)->GetName()));
+      paveText->AddText("|#eta| < 0.9");
+      
+      gPad->SetLeftMargin(0.15);
+      hist->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
+      hist->GetXaxis()->SetTitleOffset(1.5);
+      hist->GetYaxis()->SetTitleOffset(1.7);
+      hist->SetStats(0);
+      hist->SetTitle("a) Correlation");
+      TH2* clone = (TH2*) hist->Clone(Form("%s_clone", hist->GetName()));
+    //   clone->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
+      clone->Draw("SURF1");
+      paveText->Draw();
+    }
+  }
+}
+
 Bool_t disableUncertainties = kTRUE;
 
-void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1** projPhi, TH1** projEta, TH1** projPhiSyst, TH1** projEtaSyst, Bool_t Ratio = kFALSE)
+void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1** projPhi = 0, TH1** projEta = 0, TH1** projPhiSyst = 0 , TH1** projEtaSyst = 0, Bool_t Ratio = kFALSE)
 {
   Float_t etaLimit = kEtaLimit;
   Float_t outerLimit = kOuterLimit;
@@ -3072,7 +3249,7 @@ void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1**
   hist->GetZaxis()->SetTitle(kCorrFuncTitle);
   hist->GetZaxis()->SetTitleOffset(1.9);
   
-  if (graphID >= 15)
+  if (graphID >= 20)
   {
     etaLimit = 0.5;
     outerLimit = 0.99;
@@ -3103,7 +3280,7 @@ void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1**
   hist->SetStats(0);
   hist->SetTitle("a) Correlation");
   TH2* clone = (TH2*) hist->Clone(Form("%s_clone", hist->GetName()));
-  clone->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
+//   clone->GetXaxis()->SetRangeUser(-TMath::Pi() / 2, TMath::Pi() / 2);
   clone->Draw("SURF1");
   paveText->Draw();
   
@@ -3180,6 +3357,9 @@ void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1**
       etaRatio->SetBinContent(x,proj2->GetBinContent(proj2->FindBin(proj2->GetBinCenter(x)))+proj2->GetBinContent(proj2->FindBin(-1*proj2->GetBinCenter(x))));
       etaRatio->SetBinError(x,sqrt((proj2->GetBinError(proj2->FindBin(proj2->GetBinCenter(x)))*proj2->GetBinError(proj2->FindBin(proj2->GetBinCenter(x)))+proj2->GetBinError(proj2->FindBin(-1*proj2->GetBinCenter(x)))*proj2->GetBinError(proj2->FindBin(-1*proj2->GetBinCenter(x))))));
     }
+    
+  if (!projPhi || ! projEta)
+    return;
 
   if (Ratio)
   {
@@ -3197,7 +3377,7 @@ void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1**
   c->SaveAs(Form("ex/%s.png", c->GetName()));
   c->SaveAs(Form("ex/%s.eps", c->GetName()));
 
-  if (1)
+  if (0)
   {
     c = new TCanvas(Form("ex_%d_%d_%d_a", i, j, histId), Form("ex_%d_%d_%d_a", i, j, histId), 800, 800);
     gPad->SetLeftMargin(0.15);
@@ -3523,7 +3703,7 @@ void DrawFullCentralityDependence(const char* histFileName)
   for (Int_t histId = 0; histId<6; histId++)
   {
     projectionsPhi[count] = 0;
-    DrawExample(histFileName, 0, 1, histId, &projectionsPhi[count], &projectionsEta[count], &projectionsPhiSyst[count], &projectionsEtaSyst[count]);
+    DrawExample(histFileName, 1, 1, histId, &projectionsPhi[count], &projectionsEta[count], &projectionsPhiSyst[count], &projectionsEtaSyst[count]);
     count++;
   }
   
@@ -4234,8 +4414,8 @@ void v2DependencyToy2(const char* DataFileName, const char* ToyFileName, const c
 
 void DrawResults()
 {
-  const char* mainFile = "graphs_120511.root";
-  const char* wingFile = "graphs_120511_wingremoved.root";
+  const char* mainFile = "graphs_131112.root";
+  const char* wingFile = "graphs_131112_wingremoved.root";
   const char* hijing = "graphs_hijing_120515.root";
   const char* ampt = "graph_ampt_120516.root";
   
@@ -4283,29 +4463,45 @@ void DrawResults()
   }
 }
 
-void DrawAMPTComparison()
+void DrawAllMCComparisons()
 {
-  const char* mainFile = "graphs_120511.root";
-  const char* wingFile = "graphs_120511_wingremoved.root";
-  const char* ampt2 = "ampt_stringmelting_norescattering/graphs.root";
-  const char* ampt3 = "ampt_default/graphs.root";
-  
-  drawLogo = 1;
+  const char* mainFile = "graphs_131112.root";
+  const char* wingFile = "graphs_131112_wingremoved.root";
+  const char* hijing = "graphs_hijing_131122.root";
+  const char* ampt2 = "graphs_ampt_norescattering_131115.root";
+  const char* ampt3 = "graphs_ampt_default_131115.root";
+  const char* ampt4 = "graphs_ampt_nostringmelting_131115.root";
+  const char* ampt5 = "graphs_ampt_tuned_131115.root";
+
+//   drawLogo = 1;
 //   DrawResultsCentrality(mainFile, wingFile); return;
 //   MCLabel = "Data  HIJING 1.36 / Pythia 6 P-0"; MCComparison(mainFile, wingFile, hijing, 0); return;
 
   gROOT->SetBatch(kTRUE);
 
-  skipGraphList[0] = 5;
-  skipGraphList[1] = 10;
+//   skipGraphList[0] = 5;
+//   skipGraphList[1] = 10;
+
+  MCLabel = "Data  HIJING / Pythia 6 P-0";
+  fgFolder.Form("results/%s", "hijing");
+  MCComparison(mainFile, wingFile, hijing, 0);
 
-  MCLabel = "Data  AMPT 2.25 no resc / Pythia 6 P-0";
+  MCLabel = "Data  AMPT melting / Pythia 6 P-0";
   fgFolder.Form("results/%s", "ampt-norescattering");
   MCComparison(mainFile, wingFile, ampt2, 0);
 
-  MCLabel = "Data  AMPT 1.25 / Pythia 6 P-0";
+  MCLabel = "Data  AMPT melting + resc / Pythia 6 P-0";
   fgFolder.Form("results/%s", "ampt-default");
   MCComparison(mainFile, wingFile, ampt3, 0);
+
+  MCLabel = "Data  AMPT resc / Pythia 6 P-0";
+  fgFolder.Form("results/%s", "ampt-nostringmelting");
+  MCComparison(mainFile, wingFile, ampt4, 0);
+
+  MCLabel = "Data  AMPT melting + resc 50% flow / Pythia 6 P-0";
+  fgFolder.Form("results/%s", "ampt-tuned");
+  MCComparison(mainFile, wingFile, ampt5, 0);
+  
 }
 
 void CalculateIAA(const char* GraphFileName, const char* HistFileName, const char* outputFileName = "graphs.root")
@@ -4642,9 +4838,9 @@ void CompareHistDraw(const char* HistFileName1, const char* HistFileName2, Int_t
     }
   }
   
-  hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
-  hist2->GetYaxis()->SetRangeUser(-1.99, 1.99);
-  ratio->GetYaxis()->SetRangeUser(-1.99, 1.99);
+  hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  hist2->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  ratio->GetYaxis()->SetRangeUser(-1.79, 1.79);
   
   c->cd(1);
   hist1->Draw("surf1");
@@ -4660,6 +4856,62 @@ void CompareHistDraw(const char* HistFileName1, const char* HistFileName2, Int_t
   hist2->ProjectionY("p2", hist2->GetXaxis()->FindBin(-1.5), hist2->GetXaxis()->FindBin(1.5))->DrawCopy("SAME")->SetLineColor(2);
 }
 
+void CompareHistDraw(const char* HistFileName1, Int_t i, Int_t j, Int_t histId, Int_t i2, Int_t j2, Int_t histId2)
+{
+  TFile* input1 = TFile::Open(HistFileName1);
+
+  TH2* hist1 = (TH2*) input1->Get(Form("dphi_%d_%d_%d", i, j+1, histId));
+  if (!hist1)
+  {
+    cout << "Hist1 does not exist." << endl;
+    return;
+  }
+
+  TH2* hist2 = (TH2*) input1->Get(Form("dphi_%d_%d_%d", i2, j2+1, histId2));
+  if (!hist2)
+  {
+    cout << "Hist2 does not exist." << endl;
+    return;
+  }
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
+  
+  Printf("Integrals: %f %f", hist1->Integral(1, hist1->GetNbinsX(), hist1->GetYaxis()->FindBin(-1.19), hist1->GetYaxis()->FindBin(1.19)), hist2->Integral(1, hist1->GetNbinsX(), hist1->GetYaxis()->FindBin(-1.19), hist1->GetYaxis()->FindBin(1.19)));
+
+  if (1) //SubtractEtaGap
+  {
+    Float_t etaLimit = kEtaLimit;
+    Float_t outerLimit = kOuterLimit;
+
+    SubtractEtaGap(hist1, etaLimit, outerLimit, kFALSE);
+    SubtractEtaGap(hist2, etaLimit, outerLimit, kFALSE);
+  }
+  
+  TCanvas *c = new TCanvas(Form("c%d", canvasCount), Form("c%d", canvasCount), 600, 900);
+  c->Divide(1,3);
+  canvasCount++;
+
+  TH2* ratio = (TH2*) hist1->Clone("ratio");
+  ratio->Divide(hist2);
+  
+  if (hist1->GetEntries() < 1e4 || hist2->GetEntries() < 1e4)
+  {
+    cout << "Hist1: " << hist1->GetEntries() << " entries" << endl << "Hist2: " << hist2->GetEntries() << " entries" << endl;
+  }
+
+  hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  hist2->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  ratio->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  
+  c->cd(1);
+  hist1->Draw("surf1");
+  c->cd(2);
+  hist2->Draw("surf1");
+  c->cd(3);
+  ratio->Draw("colz");
+}
+
 void CompareSTARpTa(const char* STARFileName, const char* GraphFileName)
 {
   TCanvas *c1 = new TCanvas("c1", "STAR comparison in dphi, pTa", 800, 600);
@@ -5550,7 +5802,7 @@ void TwoPlusOneCorrelations()
   TH2D* as2 = new TH2D("as2", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   TH2D* ns2Same = new TH2D("ns2Same", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   TH2D* as2Same = new TH2D("as2Same", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
-  TH2D* as2AllSame = new TH2D("as2AllSame", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
+//   TH2D* as2AllSame = new TH2D("as2AllSame", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   TH2D* ns2Diff = new TH2D("ns2Diff", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   TH2D* as2Diff = new TH2D("as2Diff", ";delta eta; delta phi", bins, -2, 2, bins, -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
 
@@ -5643,7 +5895,7 @@ void TwoPlusOneCorrelations()
 
        Float_t weight = 1;
        // apply flow
-//     weight *= weights[j] * weights[k];
+       weight *= weights[j] * weights[k];
        
        ns->Fill(deltaEta, deltaPhi, weight);
       }
index 6b3ad86..666b532 100644 (file)
@@ -146,7 +146,8 @@ void DivideGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
     }
 
     Float_t value = graph1->GetY()[bin1] / graph2Extrapolated;
-    Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
+//     Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
+    Float_t error = value * TMath::Sqrt(TMath::Abs(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) - TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2)));
 
     graph1->GetY()[bin1] = value;
     graph1->GetEY()[bin1] = error;
@@ -2890,7 +2891,7 @@ void DrawAsFunctionOfPt(const char* graphFile1, Int_t graphID, Int_t centrID, In
   graph->Draw("AP");
 }
 
-TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** graphFiles, Int_t graphID, Int_t centrID, Int_t* maxpTIndex, const char** titles = 0)
+TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** graphFiles, Int_t graphID, Int_t centrID, Int_t* maxpTIndex, const char** titles = 0, Int_t* centrIndex = 0)
 {
   c = new TCanvas;
   if (n > 4)
@@ -2902,7 +2903,7 @@ TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** gra
   
   for (Int_t i=0; i<n; i++)
   {
-    graph = GetGraphAsFunctionOfPt(TString(baseFile) + graphFiles[i], graphID, centrID, maxpTIndex[i]);
+    graph = GetGraphAsFunctionOfPt(TString(baseFile) + graphFiles[i], graphID, (centrIndex) ? centrIndex[i] : centrID, maxpTIndex[i]);
     graph->SetLineColor((i%4)+1);
     graph->SetMarkerColor((i%4)+1);
     graph->SetMarkerStyle(20+i);
@@ -2919,6 +2920,7 @@ TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** gra
 //     graph->Print();
     
     c->cd();
+    
     if (1 || i < 4)
       graph->Draw((i == 0) ? "AP" : "PSAME");
     else
@@ -2935,12 +2937,14 @@ TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** gra
     
     if (i >= 4)
     {
-      graph = (TGraphErrors*) graph->Clone();
-//       SetYError(graph, 0);
-      DivideGraphs(graph, allGraphs[i%4]);
+      graphclone = (TGraphErrors*) graph->Clone();
+//       SetYError(graphclone, 0);
+      DivideGraphs(graphclone, allGraphs[i%4]);
       c2->cd();
-      graph->Draw((i == 4) ? "AP" : "PSAME");
-      graph->GetYaxis()->SetRangeUser(0.5, 1.5);
+      graphclone->Draw((i == 4) ? "AP" : "PSAME");
+      graphclone->GetYaxis()->SetRangeUser(0.5, 1.5);
+      
+      GraphShiftX(graph, 0.02);
     }
   }
   
@@ -2952,7 +2956,7 @@ TGraphErrors* DrawAsFunctionOfPt(const char* baseFile, Int_t n, const char** gra
 
 void DrawAsFunctionOfPt()
 {
-  const char* baseFile = "graphs_130513b"; Float_t maxY = 0.25;
+  const char* baseFile = "graphs_130606"; Float_t maxY = 0.25;
 //   const char* baseFile = "graphs_130503_nosub"; Float_t maxY = 0.4;
   const char* graphFiles[] = { ".root", "_pions.root", "_kaons.root", "_protons.root", 
   "_wingremoved.root", "_wingremoved_pions.root", "_wingremoved_kaons.root", "_wingremoved_protons.root",
@@ -2964,22 +2968,20 @@ void DrawAsFunctionOfPt()
       "h-h wingremoved", "h-#pi wingremoved", "h-K wingremoved", "h-p wingremoved", 
       "h", "#pi", "K", "p", 
       "h wingremoved", "#pi wingremoved", "K wingremoved", "p wingremoved" };
-  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, 7, 7, 7, 7 };
+  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
   
   Int_t graphID = 4;
-  Int_t centrID = 1;
+  Int_t centrID = 0;
   
   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, titles);
   first->GetYaxis()->SetRangeUser(0, maxY);
   
-//   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, titles+4);
-//   first->GetYaxis()->SetRangeUser(0, maxY);
+  first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, titles+4);
+  first->GetYaxis()->SetRangeUser(0, maxY);
 
   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+8, graphID, centrID, max, titles+8);
   first->GetYaxis()->SetRangeUser(0, maxY);
   
-  return;
-
   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+12, graphID, centrID, max, titles+12);
   first->GetYaxis()->SetRangeUser(0, maxY);
 
@@ -3008,10 +3010,10 @@ void DrawAsFunctionOfPt2()
 void DrawAsFunctionOfPt3()
 {
 //   const char* baseFile = "graphs_130515"; Float_t maxY = 0.25;
-  const char* baseFile = "graphs_130515_nosub"; Float_t maxY = 0.4;
+  const char* baseFile = "graphs_130531_mc_nosub"; Float_t maxY = 0.4;
   const char* graphFiles[] = { ".root", "_pions.root", "_kaons.root", "_protons.root", "_unfolded.root", "_pions_unfolded.root", "_kaons_unfolded.root", "_protons_unfolded.root", "_allpt_unfolded.root", "_allpt_pions_unfolded.root", "_allpt_kaons_unfolded.root", "_allpt_protons_unfolded.root" };
   const char* titles[] = { "h-h", "h-#pi", "h-K", "h-p", "h", "#pi", "K", "p", "h allpt", "#pi allpt", "K allpt", "p allpt" };
-  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, 7, 7, 7, 7 };
+  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, 6, 6, 6, 6 };
   
   Int_t graphID = 4;
   
@@ -3024,7 +3026,7 @@ void DrawAsFunctionOfPt3()
   first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, 0, max, titles);
   first->GetYaxis()->SetRangeUser(0, maxY);
   
-  return;
+//   return;
 
   first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+8, graphID, 0, max+8, titles+8);
   first->GetYaxis()->SetRangeUser(0, maxY);
@@ -3036,10 +3038,10 @@ void DrawAsFunctionOfPt3()
 void DrawAsFunctionOfPt4()
 {
 //   const char* baseFile = "graphs_130503";
-  const char* baseFile = "graphs_130503_nosub";
-  const char* graphFiles[] = { "_allpt_unfolded.root", "_allpt_pions_unfolded.root", "_allpt_kaons_unfolded.root", "_allpt_protons_unfolded.root", "_allpt.root", "_allpt_pions.root", "_allpt_kaons.root", "_allpt_protons.root" };
+  const char* baseFile = "graphs_130618";
+  const char* graphFiles[] = { "_unfolded.root", "_pions_unfolded.root", "_kaons_unfolded.root", "_protons_unfolded.root", "_otherbaseline_unfolded.root", "_otherbaseline_pions_unfolded.root", "_otherbaseline_kaons_unfolded.root", "_otherbaseline_protons_unfolded.root" };
 
-  Int_t max[8] = { 8, 8, 8, 8, 8, 8, 8, 8 };
+  Int_t max[8] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
   
   Int_t centrID = 0;
   
@@ -3076,6 +3078,29 @@ void DrawAsFunctionOfPt5()
   first->GetYaxis()->SetRangeUser(0, maxY); 
 }
 
+void DrawAsFunctionOfPt7()
+{
+  const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130604"; Float_t maxY = 0.4;
+  const char* graphFiles[] = { 
+    "_Hadrons.root", "_Pions.root", "_Kaons.root", "_Protons.root",
+    "_Hadrons_NoEff.root", "_Pions_NoEff.root", "_Kaons_NoEff.root", "_Protons_NoEff.root"
+  };
+  
+  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
+  
+  Int_t graphID = 4;
+  Int_t centrID = 0;
+  
+  first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, centrID, max, 0);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+
+  first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, centrID, max, 0);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+
+  first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, centrID, max, 0);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+}
+
 void DrawAsFunctionOfPt6()
 {
   const char* baseFile = "graphs_130515b_nosub"; Float_t maxY = 0.4;
@@ -3099,6 +3124,29 @@ void DrawAsFunctionOfPt6()
   first->GetYaxis()->SetRangeUser(0, maxY); 
 }
 
+void CompareCentralities()
+{
+  const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130604"; Float_t maxY = 0.4;
+  const char* graphFiles[] = { 
+    "_Hadrons.root", "_Pions.root", "_Kaons.root", "_Protons.root",
+    "_Hadrons.root", "_Pions.root", "_Kaons.root", "_Protons.root"
+  };
+  
+  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
+  Int_t centrID[] = { 0, 0, 0, 0, 1, 1, 1, 1 };
+  
+  Int_t graphID = 4;
+  
+  first = DrawAsFunctionOfPt(baseFile, 4, graphFiles, graphID, 0, max, 0, centrID);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+
+  first = DrawAsFunctionOfPt(baseFile, 4, graphFiles+4, graphID, 0, max, 0, centrID+4);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+
+  first = DrawAsFunctionOfPt(baseFile, 8, graphFiles, graphID, 0, max, 0, centrID);
+  first->GetYaxis()->SetRangeUser(0, maxY); 
+}
+
 void CompareCentralityEstimators()
 {
   const char* baseFile = "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_"; Float_t maxY = 0.25;
@@ -3153,6 +3201,24 @@ void DrawLikeUnlikeSign()
   first->GetYaxis()->SetRangeUser(0, maxY);
 }
 
+void Drawv1()
+{
+  Float_t maxY = 0.35;
+  const char* graphFiles[] = { 
+    "graphs_fwdA_130801.root", "graphs_fwdC_130801.root", "~/Dropbox/pid-ridge/pPb/graphs_LHC13bc_20130709_Hadrons.root"
+  };
+  
+   const char* titles[] = { 
+      "A", "C", "CB" };
+  Int_t max[] = { NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs, NGraphs };
+  
+  Int_t graphID = 12;
+  Int_t centrID = 0;
+  
+  first = DrawAsFunctionOfPt("", 3, graphFiles, graphID, centrID, max, titles);
+  first->GetYaxis()->SetRangeUser(0, maxY);
+}
+
 void DrawRelativeError(Int_t id = 0)
 {
   if (id == 0)
@@ -3221,7 +3287,7 @@ void DrawRelativeError(Int_t id = 0)
   graph->DrawClone("PSAME");
   legend->AddEntry(graph->DrawClone("PSAME"), "comb", "P");
   
-  line = new TLine(0, 0.1, 4.0, 0.1);
+  line = new TLine(0, 0.05, 4.0, 0.05);
   line->Draw();
   legend->AddEntry(line, "K0-Kch unc", "L");
   legend->SetHeader(label);  
index 94dcb09..a8d9699 100644 (file)
@@ -1,5 +1,5 @@
-// 0       1       2       3       4  5  6             7     8             9                            10         11
-// nsyield,asyield,nswidth,aswidth,v2,v3,nsasyieldratio,v3/v2,remainingpeak,remainingjetyield/ridgeyield,chi2(v2v3),baseline
+// 0       1       2       3       4  5  6             7     8             9                            10         11       12
+// nsyield,asyield,nswidth,aswidth,v2,v3,nsasyieldratio,v3/v2,remainingpeak,remainingjetyield/ridgeyield,chi2(v2v3),baseline,v1
 const char* graphTitles[] = { "NS Ridge Yield", "AS Ridge Yield", "NS Width", "AS Width", "v2", "v3", "NS Yield / AS Yield", "v3 / v2", "remaining peak in %", "remaining jet / NS yield in %", "chi2/ndf", "baseline", "v1" };
 const Int_t NGraphs = 10; // pt index  //10 for the paper
 const Int_t NHists = 13; 
@@ -983,16 +983,16 @@ void CorrelationSubtractionHistogram(const char* fileName, Bool_t centralityAxis
     Int_t n = 10; //fine binning with low pt TPC only
     Int_t is[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
     Int_t js[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
-    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
+    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
   }
   
   if (gBinning == 1)
   {
     // for extraction when trigger pT is from whole range
-    Int_t n = 8;
-    Int_t is[] =    { 0, 0, 0, 0, 0, 0, 0, 0 };
-    Int_t js[] =    { 1, 2, 3, 4, 5, 6, 7, 9 };
-    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1, 1 };
+    Int_t n = 6;
+    Int_t is[] =    { 0, 0, 0, 0, 0, 0, 0 };
+    Int_t js[] =    { 1, 2, 3, 4, 5, 6, 9 };
+    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1 };
   }
 
   if (gBinning == 2)
@@ -1017,7 +1017,7 @@ void CorrelationSubtractionHistogram(const char* fileName, Bool_t centralityAxis
   if (gBinning == 4)
   {
     // default with all but only symmetric bins
-    Int_t n = 7;
+    Int_t n = 6;
     // sorted by pT,T
     Int_t is[] =    { 0, 1, 2, 3, 4, 5, 8 };
     Int_t js[] =    { 1, 2, 3, 4, 5, 6, 9 };
@@ -1131,12 +1131,13 @@ void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr,
   TH2* hist2 = (TH2*) fin->Get(Form("dphi_%d_%d_%d", i, j, ((gStudySystematic == 60) ? 6 : 4)));
   TH1* refMult = (TH1*) fin->Get("refMult");
   
+  if (!hist1 || !hist2)
+    return;
+
   hist1 = (TH2*) hist1->Clone();
   hist2 = (TH2*) hist2->Clone();
   refMult = (TH1*) refMult->Clone();
   
-  if (!hist1 || !hist2)
-    return 0;
   // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
   hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
   hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
@@ -1587,6 +1588,21 @@ void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr,
   Float_t baseLine2 = hist2->IntegralAndError(phi1Z, phi2Z, eta1, eta6, baseLineE2);
   baseLine2 /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
   baseLineE2 /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
+
+  // get baseline of peripheral bin with large eta gap
+  Int_t eta1Alternative = hist1->GetYaxis()->FindBin(-etaMax + 0.001);
+  Int_t eta2Alternative = hist1->GetYaxis()->FindBin(-1.2 - 0.001);
+  Int_t eta3Alternative = hist1->GetYaxis()->FindBin(1.2 + 0.001);
+  Int_t eta6Alternative = hist1->GetYaxis()->FindBin(etaMax - 0.001);
+  Int_t phi1ZAlternative = hist1->GetXaxis()->FindBin(-TMath::Pi()/2 + 0.001);
+  Int_t phi2ZAlternative = hist1->GetXaxis()->FindBin( TMath::Pi()/2 - 0.001);
+  
+  Float_t baseLine2NSGapE, baseLine2NSGapE2;
+  Float_t baseLine2NSGap = hist2->IntegralAndError(phi1ZAlternative, phi2ZAlternative, eta1Alternative, eta2Alternative, baseLine2NSGapE);
+  baseLine2NSGap += hist2->IntegralAndError(phi1ZAlternative, phi2ZAlternative, eta3Alternative, eta6Alternative, baseLine2NSGapE2);
+  baseLine2NSGapE = TMath::Sqrt(baseLine2NSGapE * baseLine2NSGapE + baseLine2NSGapE2 * baseLine2NSGapE2);
+  baseLine2NSGap /= (phi2ZAlternative - phi1ZAlternative + 1) * (eta2Alternative - eta1Alternative + 1 + eta6Alternative - eta3Alternative + 1);
+  baseLine2NSGapE /= (phi2ZAlternative - phi1ZAlternative + 1) * (eta2Alternative - eta1Alternative + 1 + eta6Alternative - eta3Alternative + 1);
   
   if(gDoSubtraction)hist1->Add(hist2, -1);
   
@@ -2037,12 +2053,22 @@ void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr,
     if (v2v3->GetParameter(1) > 0)
     {
       if(gDoSubtraction){
-       v2value = TMath::Sqrt(v2v3->GetParameter(1) / (baseLine + diffMinParam0));
-       v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + baseLineE * baseLineE / baseLine / baseLine);
+       if (1)
+       {
+         v2value = TMath::Sqrt(v2v3->GetParameter(1) / (baseLine + diffMinParam0));
+         v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + baseLineE * baseLineE / baseLine / baseLine);
+         Printf("-> v2 = %f +- %f (%f %f = %f)", v2value, v2E, baseLine, diffMinParam0, baseLine + diffMinParam0);
+       }
+       else
+       {
+         v2value = TMath::Sqrt(v2v3->GetParameter(1) / (v2v3->GetParameter(0) + baseLine2NSGap));
+         v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + v2v3->GetParError(0) * v2v3->GetParError(0) / v2v3->GetParameter(0) / v2v3->GetParameter(0) + baseLine2NSGapE * baseLine2NSGapE / baseLine2NSGap / baseLine2NSGap);
+         Printf("-> v2 = %f +- %f (%f %f = %f)", v2value, v2E, v2v3->GetParameter(0), baseLine2NSGap, v2v3->GetParameter(0) + baseLine2NSGap);
+       }
       }else{
        v2value = TMath::Sqrt(v2v3->GetParameter(1) / v2v3->GetParameter(0));
        v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + v2v3->GetParError(0) * v2v3->GetParError(0) / v2v3->GetParameter(0) / v2v3->GetParameter(0));
-      }        
+      }       
     }
     if (symmetricpT)//all points filled (even if v2v3->GetParameter(1) < 0), easier to plot after (same point number for different particles)
       AddPoint(graph[4], xValue1vn, v2value, 0, v2E);
@@ -3249,7 +3275,7 @@ void DivideOutReferenceParticle(const char* graphFile, const char* outputFile =
 
   for (Int_t graphID = 4; graphID <= 4; graphID++) //v2 + v3
   {
-    const Int_t posFullPt = 7; // position where vn for unfolding is located (i.e. full pT,a and pT,t range)
+    const Int_t posFullPt = 6; // position where vn for unfolding is located (i.e. full pT,a and pT,t range)
     Bool_t allPt = TString(graphFile).Contains("allpt");
     Int_t maxPt = (allPt) ? posFullPt : NGraphs;
     for (Int_t i=0; i<maxPt; i++)
@@ -3308,15 +3334,15 @@ void ExtractForSkalarProductComparisonAllSpecies()
 {
   gROOT->SetBatch(kTRUE);
 
-  TString baseInput("dphi_corr_130515");
-  TString baseOutput("graphs_130515b");
-  baseOutput += "_nosub";
+  TString baseInput("dphi_corr_LHC13bc_20130604");
+  TString baseOutput("graphs_130618_otherbaseline");
+//   baseOutput += "_nosub";
 //   TString postfix("_wingremoved");
   TString postfix;
   
   if (1)
   {
-    gBinning = 4;
+    gBinning = 0;
     gReferenceFile = baseOutput + ".root";
     ExtractForSkalarProductComparison(baseInput + "_Hadrons" + postfix, baseOutput, kTRUE); 
     ExtractForSkalarProductComparison(baseInput + "_Pions" + postfix, baseOutput + "_pions", kTRUE);
@@ -3327,14 +3353,13 @@ void ExtractForSkalarProductComparisonAllSpecies()
   if (0)
   {
     gBinning = 1;
-    baseInput = "dphi_corr_130506";
     baseInput += "_allpt";
     baseOutput += "_allpt";
     gReferenceFile = baseOutput + ".root";
-    ExtractForSkalarProductComparison(baseInput + "_Hadrons" + postfix, baseOutput, kTRUE); 
-    ExtractForSkalarProductComparison(baseInput + "_Pions" + postfix, baseOutput + "_pions", kTRUE);
-    ExtractForSkalarProductComparison(baseInput + "_Protons" + postfix, baseOutput + "_protons", kTRUE);
-    ExtractForSkalarProductComparison(baseInput + "_Kaons" + postfix, baseOutput + "_kaons", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_55_Hadrons_symmetrized" + postfix, baseOutput, kTRUE); 
+    ExtractForSkalarProductComparison(baseInput + "_63_Pions" + postfix, baseOutput + "_pions", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_63_Protons" + postfix, baseOutput + "_protons", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_63_Kaons" + postfix, baseOutput + "_kaons", kTRUE);
   }
 }
 
diff --git a/PWGCF/Correlations/macros/dphicorrelations/pA_PID_Fwd.C b/PWGCF/Correlations/macros/dphicorrelations/pA_PID_Fwd.C
new file mode 100644 (file)
index 0000000..f5d4a6c
--- /dev/null
@@ -0,0 +1,3006 @@
+// 0       1       2       3       4  5  6             7     8             9                            10         11
+// nsyield,asyield,nswidth,aswidth,v2,v3,nsasyieldratio,v3/v2,remainingpeak,remainingjetyield/ridgeyield,chi2(v2v3),baseline
+const char* graphTitles[] = { "NS Ridge Yield", "AS Ridge Yield", "NS Width", "AS Width", "v2", "v3", "NS Yield / AS Yield", "v3 / v2", "remaining peak in %", "remaining jet / NS yield in %", "chi2/ndf", "baseline", "v1" };
+const Int_t NGraphs = 28; // pt index
+const Int_t NHists = 13; 
+TGraphErrors*** graphs = 0;
+
+const char* kCorrFuncTitle = "#frac{1}{#it{N}_{trig}} #frac{d^{2}#it{N}_{assoc}}{d#Delta#etad#Delta#varphi} (rad^{-1})";
+// const char* kCorrFuncTitle = "1/N_{trig} dN_{assoc}/d#Delta#etad#Delta#varphi (1/rad)";
+const char* kProjYieldTitlePhi = "1/#it{N}_{trig} d#it{N}_{assoc}/d#Delta#varphi per #Delta#eta (rad^{-1})";
+const char* kProjYieldTitleEta = "1/#it{N}_{trig} d#it{N}_{assoc}/d#Delta#eta per #Delta#varphi (rad^{-1})";
+// const char* kProjYieldTitlePhiOrEta = "1/N_{trig} dN_{assoc}/d#Delta#varphi (1/rad) , dN_{assoc}/d#Delta#eta";
+const char* kProjYieldTitlePhiOrEta = "#frac{1}{#it{N}_{trig}} #frac{d#it{N}_{assoc}}{d#Delta#varphi} (rad^{-1}) , d#it{N}_{assoc}/d#Delta#eta";
+
+Float_t etaMax = 4.0;
+Bool_t gUseAnalyticalVn=0;
+Bool_t gDoSubtraction=1;
+Bool_t gDoEtaGap=1;
+Bool_t gDoEtaGapAS=0;
+Bool_t gExcludeNearSidePeakFromPhiProj=0;
+Int_t gBinning = 0;
+
+void CreateGraphStructure()
+{
+  graphs = new TGraphErrors**[NGraphs];
+  for (Int_t i=0; i<NGraphs; i++)
+    {
+      graphs[i] = new TGraphErrors*[NHists];
+      for (Int_t j=0; j<NHists; j++)
+       {
+         graphs[i][j] = new TGraphErrors;
+         graphs[i][j]->GetYaxis()->SetTitle(graphTitles[j]);
+       }
+    }
+}
+
+void WriteGraphs(const char* outputFileName = "graphs.root")
+{
+  TFile::Open(outputFileName, "RECREATE");
+  for (Int_t i=0; i<NGraphs; i++)
+    for (Int_t j=0; j<NHists; j++)
+      {
+       if (!graphs[i][j])
+         continue;
+       graphs[i][j]->GetYaxis()->SetTitle(graphTitles[j]);
+       graphs[i][j]->Write(Form("graph_%d_%d", i, j));
+      }
+
+  gFile->Close();
+}
+
+void ReadGraphs(const char* fileName = "graphs.root")
+{
+  CreateGraphStructure();
+  TFile* file = TFile::Open(fileName);
+  for (Int_t i=0; i<NGraphs; i++)
+    for (Int_t j=0; j<NHists; j++)
+      graphs[i][j] = (TGraphErrors*) file->Get(Form("graph_%d_%d", i, j));
+}
+
+void AddPoint(TGraphErrors* graph, Float_t x, Float_t y, Float_t xe, Float_t ye)
+{
+  graph->SetPoint(graph->GetN(), x, y);
+  graph->SetPointError(graph->GetN() - 1, xe, ye);
+}
+
+void DrawLatex(Float_t x, Float_t y, Int_t color, const char* text, Float_t textSize = 0.06)
+{
+  TLatex* latex = new TLatex(x, y, text);
+  latex->SetNDC();
+  latex->SetTextSize(textSize);
+  latex->SetTextColor(color);
+  latex->Draw();
+}
+
+void PadFor2DCorr()
+{
+  gPad->SetPad(0, 0, 1, 1);
+  gPad->SetLeftMargin(0.2);
+  gPad->SetTopMargin(0.05);
+  gPad->SetRightMargin(0.05);
+  gPad->SetBottomMargin(0.05);
+  gPad->SetTheta(61.62587);
+  gPad->SetPhi(45);
+}
+
+void DivideGraphs(TGraphErrors* graph1, TGraphErrors* graph2)
+{
+  //   graph1->Print();
+  //   graph2->Print();
+
+  graph1->Sort();
+  graph2->Sort();
+
+  for (Int_t bin1 = 0; bin1 < graph1->GetN(); bin1++)
+    {
+      Float_t x = graph1->GetX()[bin1];
+
+      Int_t bin2 = 0;
+      for (bin2 = 0; bin2<graph2->GetN(); bin2++)
+       if (graph2->GetX()[bin2] >= x)
+         break;
+
+      if (bin2 == graph2->GetN())
+       bin2--;
+
+      if (bin2 > 0)
+       if (TMath::Abs(graph2->GetX()[bin2-1] - x) < TMath::Abs(graph2->GetX()[bin2] - x))
+         bin2--;
+
+      if (graph1->GetY()[bin1] == 0 || graph2->GetY()[bin2] == 0 || bin2 == graph2->GetN())
+       {
+         Printf("%d %d removed", bin1, bin2);
+         graph1->RemovePoint(bin1--);
+         continue;
+       }
+
+      Float_t graph2Extrapolated = graph2->GetY()[bin2];
+      if (TMath::Abs(x - graph2->GetX()[bin2]) > 0.0001)
+       {
+         Printf("%f %f %d %d not found", x, graph2->GetX()[bin2], bin1, bin2);
+         graph1->RemovePoint(bin1--);
+         continue;
+       }
+
+      Float_t value = graph1->GetY()[bin1] / graph2Extrapolated;
+      Float_t error = value * TMath::Sqrt(TMath::Power(graph1->GetEY()[bin1] / graph1->GetY()[bin1], 2) + TMath::Power(graph2->GetEY()[bin2] / graph2->GetY()[bin2], 2));
+
+      graph1->GetY()[bin1] = value;
+      graph1->GetEY()[bin1] = error;
+
+      //     Printf("%d %d %f %f %f %f", bin1, bin2, x, graph2Extrapolated, value, error);
+    }
+}
+
+void GraphShiftX(TGraphErrors* graph, Float_t offset)
+{
+  for (Int_t i=0; i<graph->GetN(); i++)
+    graph->GetX()[i] += offset;
+}
+
+TGraphErrors* ReadHepdata(const char* fileName, Bool_t errorsAreAdded = kFALSE, Int_t skipYErrors = 0, Int_t skipXerrors = 1)
+{
+  // expected format: x [x2] y [ye] [ye2] [xe]
+  //
+  // skipYErrors:   0 --> ye present
+  //                1 --> no errors ye
+  //                2 --> y and ye are lower and upper error, i.e. y' = (y + ye) / 2 and ye = (ye - y) / 2
+  //                3 --> ye and ye2 are stat and syst error, will be added in quadrature
+  // 
+  // skipXerrors:   0 --> xe present
+  //                1 --> no errors xe
+  //                2 --> x2 present, xe not present and is calculated from x2 - x
+  
+  ifstream fin(fileName);
+
+  graph = new TGraphErrors(0);
+
+  Double_t sum = 0;
+
+  while (fin.good())
+    {
+      char buffer[2000];
+      if (fin.peek() == '#')
+       {
+         fin.getline(buffer, 2000);
+         continue;
+       }
+  
+      Double_t x = -1;
+      Double_t x2 = -1;
+      Double_t y = -1;
+      Double_t ye = 0;
+      Double_t xe = 0;
+
+      fin >> x;
+    
+      if (skipXerrors == 2)
+       {
+         fin >> x2;
+         xe = (x2 - x + 1) / 2;
+         x = x + (x2 - x) / 2;
+       }
+    
+      fin >> y;
+
+      if (y == -1)
+       continue;
+
+      if (skipYErrors == 0)
+       {
+         ye = -1;
+         fin >> ye;
+         if (ye == -1)
+           continue;
+       }
+      else if (skipYErrors == 2)
+       {
+         ye = -1;
+         fin >> ye;
+         if (ye == -1)
+           continue;
+      
+         Double_t newy = (y + ye) / 2;
+         ye = (ye - y) / 2;
+         y = newy;
+       }
+      else if (skipYErrors == 3)
+       {
+         ye = -1;
+         fin >> ye;
+         if (ye == -1)
+           continue;
+      
+         Double_t ye2 = -1;
+         fin >> ye2;
+         if (ye2 == -1)
+           continue;
+
+         ye = TMath::Sqrt(ye*ye + ye2*ye2);
+       }
+
+      if (skipXerrors == 0)
+       {
+         xe = -1;
+         fin >> xe;
+         if (xe == -1)
+           continue;
+       }
+
+      //Printf("%f %f %f %f", x, y, xe, ye);
+
+      if (errorsAreAdded)
+       ye -= y;
+
+      graph->SetPoint(graph->GetN(), x, y);
+      graph->SetPointError(graph->GetN()-1, xe, ye);
+
+      sum += y;
+    
+      // read rest until end of line...
+      fin.getline(buffer, 2000);
+    }
+  fin.close();
+
+  Printf("%s: %f", fileName, sum);
+
+  return graph;
+}
+
+TH2* SubtractEtaGapNS(TH2* hist, Float_t etaLimit, Float_t outerLimit, Bool_t drawEtaGapDist = kFALSE)
+{
+  TString histName(hist->GetName());
+  Int_t etaBins = 0;
+
+  TH1D* etaGap = hist->ProjectionX(histName + "_1", TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)), hist->GetYaxis()->FindBin(-etaLimit - 0.01));
+  //   Printf("%f", etaGap->GetEntries());
+  if (etaGap->GetEntries() > 0)
+    etaBins += hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
+
+  TH1D* tracksTmp = hist->ProjectionX(histName + "_2", hist->GetYaxis()->FindBin(etaLimit + 0.01), TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)));
+  //   Printf("%f", tracksTmp->GetEntries());
+  if (tracksTmp->GetEntries() > 0)
+    etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
+  
+  etaGap->Add(tracksTmp);
+
+  // get per bin result
+  if (etaBins > 0)
+    etaGap->Scale(1.0 / etaBins);
+  for (Int_t i=1; i<=etaGap->GetNbinsX()/2; i++)
+    {
+      //     Printf("%d -> %d", i, etaGap->GetNbinsX()+1-i);
+      etaGap->SetBinContent(etaGap->GetNbinsX()+1-i, etaGap->GetBinContent(i));
+      etaGap->SetBinError(etaGap->GetNbinsX()+1-i, etaGap->GetBinError(i));
+    }
+  
+  if (drawEtaGapDist)
+    {
+      TH1D* centralRegion = hist->ProjectionX(histName + "_3", hist->GetYaxis()->FindBin(-etaLimit + 0.01), hist->GetYaxis()->FindBin(etaLimit - 0.01));
+    
+      //    centralRegion->Scale(1.0 / (hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
+      centralRegion->Scale(hist->GetXaxis()->GetBinWidth(1));
+
+      TCanvas* c = new TCanvas("SubtractEtaGap", "SubtractEtaGap", 800, 800);
+      gPad->SetLeftMargin(0.13);
+      centralRegion->SetStats(0);
+      TString label(centralRegion->GetTitle());
+      label.ReplaceAll(".00", " GeV/#it{c}");
+      label.ReplaceAll(".0", " GeV/#it{c}");
+      centralRegion->SetTitle(label);
+      centralRegion->SetLineColor(3);
+      centralRegion->Draw();
+      //     centralRegion->GetYaxis()->SetTitle(kProjYieldTitlePhi);
+      centralRegion->GetYaxis()->SetTitleOffset(1.6);
+      TH1* copy = etaGap->DrawCopy("SAME");
+      copy->Scale(hist->GetXaxis()->GetBinWidth(1));
+      copy->Scale((hist->GetYaxis()->FindBin(etaLimit - 0.01) - hist->GetYaxis()->FindBin(-etaLimit + 0.01) + 1));
+      copy->SetLineColor(2);
+      TLegend* legend = new TLegend(0.41, 0.73, 0.69, 0.85);
+      legend->SetFillColor(0);
+      legend->AddEntry(centralRegion, Form("|#Delta#eta| < %.1f", etaLimit), "L");
+      legend->AddEntry(copy, Form("%.1f < |#Delta#eta| < %.1f (scaled)", etaLimit, outerLimit), "L");
+      legend->Draw();
+    
+      //     DrawLatex(0.705, 0.62, 1, "Pb-Pb 2.76 TeV", 0.025);
+      //     DrawLatex(0.705, 0.58, 1, "Stat. unc. only", 0.025);
+    }
+  
+  //   new TCanvas; etaGap->DrawCopy();
+  
+  TH2* histTmp2D = (TH2*) hist->Clone("histTmp2D");
+  histTmp2D->Reset();
+  
+  for (Int_t xbin=1; xbin<=histTmp2D->GetNbinsX(); xbin++)
+    for (Int_t y=1; y<=histTmp2D->GetNbinsY(); y++)
+      histTmp2D->SetBinContent(xbin, y, etaGap->GetBinContent(xbin));
+    
+  hist->Add(histTmp2D, -1);  
+  return histTmp2D;
+}
+
+TH1* GetProjections(Int_t i, Int_t j, Int_t centr, char** label, Float_t etaBegin = 1.0, TH1** etaProj = 0)
+{
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+  if (!hist1)
+    return 0;
+  hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), etaBegin));
+
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  
+  hist1->Rebin2D(2, 1); hist1->Scale(0.5);
+  
+  //   new TCanvas; hist1->Draw("surf1");
+  
+  if (etaBegin > 0)
+    SubtractEtaGapNS(hist1, etaBegin, etaMax, kTRUE);
+  
+  tokens = TString(hist1->GetTitle()).Tokenize("-");
+  centralityStr = new TString;
+  if (tokens->GetEntries() > 2)
+    *centralityStr = tokens->At(2)->GetName();
+  if (tokens->GetEntries() > 3)
+    *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
+  *label = centralityStr->Data();
+  //   Printf("%s", label);
+  
+  proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
+  proj1x->GetXaxis()->SetTitleOffset(1);
+  proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
+
+  proj1y = ((TH2*) hist1)->ProjectionY(Form("proj1y_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
+  proj1y->Scale(hist1->GetXaxis()->GetBinWidth(1));
+  proj1y->GetXaxis()->SetTitleOffset(1);
+  proj1y->GetXaxis()->SetRangeUser(-1.99, 1.99);
+  Float_t etaPhiScale = 1.0 * (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1) / (hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1);
+  
+  if (gStudySystematic == 20)
+    {
+      Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
+      file2 = TFile::Open("non_closure.root");
+      non_closure = (TH1*) file2->Get(Form("non_closure_all_%d_%d_%d", i, j, 0));
+      for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
+       non_closure->SetBinError(bin, 0);
+    
+      proj1x->Multiply(non_closure);  
+    }
+  
+  Float_t zyam = 0;
+  if (0)
+    {  
+      clone = (TH1*) proj1x->Clone();
+      clone->Fit("pol0", "0", "", TMath::Pi()/2 - 0.2, TMath::Pi()/2);
+      zyam = clone->GetFunction("pol0")->GetParameter(0);
+    }
+  else
+    {
+      //     zyam = (proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2)) + proj1x->GetBinContent(proj1x->FindBin(-TMath::Pi()/2))) / 2;
+      //     zyam = proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2));
+      zyam = proj1x->GetBinContent(proj1x->FindBin(1.3));
+      //     zyam = proj1x->GetMinimum();
+    }
+    
+  proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
+  proj1y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
+  
+  proj1x->Scale(1.0 / (2.0 * etaMax));
+  
+  if (etaProj != 0)
+    *etaProj = proj1y;
+  return proj1x;
+}
+
+TH1* GetProjectionsNew(Int_t i, Int_t j, Int_t centr, char** label, Float_t etaBegin = 1.0, TH1** etaProj = 0)
+{
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+  if (!hist1)
+    return 0;
+  hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), etaBegin));
+
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  
+  hist1->Rebin2D(2, 1); hist1->Scale(0.5);
+  
+  //   new TCanvas; hist1->Draw("surf1");
+  
+  tokens = TString(hist1->GetTitle()).Tokenize("-");
+  centralityStr = new TString;
+  if (tokens->GetEntries() > 2)
+    *centralityStr = tokens->At(2)->GetName();
+  if (tokens->GetEntries() > 3)
+    *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
+  *label = centralityStr->Data();
+  //   Printf("%s", label);
+  
+  if (etaBegin > 0)
+    {
+      proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
+      proj1x->GetXaxis()->SetTitleOffset(1);
+      proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
+
+      proj1xR1 = ((TH2*) hist1)->ProjectionX(Form("proj2x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(-etaBegin-0.01));
+      proj1xR2 = ((TH2*) hist1)->ProjectionX(Form("proj3x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(etaBegin+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
+      proj1xR1->Add(proj1xR2);
+    
+      proj1xR1->GetXaxis()->SetTitleOffset(1);
+      proj1xR1->Scale(hist1->GetYaxis()->GetBinWidth(1));
+    
+      proj1xR1->Scale((1.0 * hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1) / (hist1->GetYaxis()->FindBin(-etaBegin-0.01) - hist1->GetYaxis()->FindBin(-etaMax+0.01) + 1 + hist1->GetYaxis()->FindBin(etaMax-0.01) - hist1->GetYaxis()->FindBin(etaBegin+0.01) + 1));
+    
+      // mirror
+      for (Int_t i=1; i<=proj1xR1->GetNbinsX()/2; i++)
+       {
+         //     Printf("%d -> %d", i, etaGap->GetNbinsX()+1-i);
+         proj1xR1->SetBinContent(proj1xR1->GetNbinsX()+1-i, proj1xR1->GetBinContent(i));
+         proj1xR1->SetBinError(proj1xR1->GetNbinsX()+1-i, proj1xR1->GetBinError(i));
+       }
+    
+      proj1x->Add(proj1xR1, -1);
+    
+      proj1x->Scale(1.0 / (2.0 * etaMax));
+  
+      //     new TCanvas; proj1xR1->DrawCopy();
+    }
+  else
+    {
+      proj1x = ((TH2*) hist1)->ProjectionX(Form("proj1x_%d_%d_%d_%.1f", i, j, centr, etaBegin), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
+      proj1x->GetXaxis()->SetTitleOffset(1);
+      proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
+
+      proj1x->Scale(1.0 / (2.0 * etaMax));
+    }
+
+  if (gStudySystematic == 20)
+    {
+      Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
+      file2 = TFile::Open("non_closure.root");
+      non_closure = (TH1*) file2->Get(Form("non_closure_all_%d_%d_%d", i, j, 0));
+      for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
+       non_closure->SetBinError(bin, 0);
+    
+      proj1x->Multiply(non_closure);  
+    }
+  
+  Float_t zyam = 0;
+  if (0)
+    {  
+      clone = (TH1*) proj1x->Clone();
+      clone->Fit("pol0", "0", "", TMath::Pi()/2 - 0.2, TMath::Pi()/2);
+      zyam = clone->GetFunction("pol0")->GetParameter(0);
+    }
+  else
+    {
+      //     zyam = (proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2)) + proj1x->GetBinContent(proj1x->FindBin(-TMath::Pi()/2))) / 2;
+      //     zyam = proj1x->GetBinContent(proj1x->FindBin(TMath::Pi()/2));
+      zyam = proj1x->GetBinContent(proj1x->FindBin(1.3));
+      //     zyam = proj1x->GetMinimum();
+    }
+    
+  proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
+  
+  return proj1x;
+}
+
+void DrawEtaDep(const char* fileName, Int_t i, Int_t j, Int_t centr)
+{
+  c = new TCanvas;
+  gPad->SetGridx();
+  gPad->SetGridy();
+  
+  TFile::Open(fileName);
+
+  legend = new TLegend(0.8, 0.6, 0.99, 0.99);
+  legend->SetFillColor(0);
+
+  for (Int_t n=0; n<4; n++)
+    {
+      Float_t eta = 0.8 + n * 0.2;
+      const char* label = 0;
+      TH1* hist = GetProjections(i, j, centr, &label, eta);
+      hist->SetStats(0);       
+      hist->SetLineColor(n+1);
+      hist->Draw((n == 0) ? "" : "SAME");
+      legend->AddEntry(hist, Form("|#Delta#eta| > %.1f", eta));
+    }
+  legend->Draw();
+}
+
+
+void DrawProjections(const char* fileName, Int_t i, Int_t j, Float_t eta = 1.0, Bool_t etaPhi = kFALSE)
+{
+  c = new TCanvas;
+  gPad->SetGridx();
+  gPad->SetGridy();
+  
+  TFile::Open(fileName);
+
+  legend = new TLegend(0.8, 0.6, 0.99, 0.99);
+  legend->SetFillColor(0);
+  
+  Int_t colors[] = { 1, 2, 1, 4, 6, 2 };
+  Int_t markers[] = { 1, 1, 5, 1, 1, 5 };
+  
+  TH1* first = 0;
+  
+  Float_t min = 100;
+  Float_t max = -100;
+  
+  for (Int_t centr=0; centr<6; centr++)
+    {
+      /*    if (centr >= 5)
+           continue;*/
+      const char* label = 0;
+      TH1* etaHist = 0;
+      TH1* hist = GetProjectionsNew(i, j, centr, &label, eta, &etaHist);
+      if (!hist)
+       continue;
+      if (etaPhi)
+       hist = etaHist;
+      hist->SetStats(0);
+      hist->SetLineColor(colors[centr]);
+      hist->SetMarkerColor(colors[centr]);
+      hist->SetMarkerStyle(markers[centr]);
+      if (etaPhi)
+       hist->GetXaxis()->SetRangeUser(-1.79, 1.79);
+      c->cd();
+      hist->Draw((centr == 0) ? "" : "SAME");
+      min = TMath::Min(min, hist->GetMinimum());
+      max = TMath::Max(max, hist->GetMaximum());
+      if (!first)
+       first = hist;
+      legend->AddEntry(hist, label);
+      //     break;
+    }
+  first->GetYaxis()->SetRangeUser(min * 1.1, max * 1.1);
+  legend->Draw();
+  
+  c->SaveAs(Form("%s_%d_%d.png", (etaPhi) ? "eta" : "phi", i, j));
+  
+  return;
+  
+  c = new TCanvas;
+  gPad->SetGridx();
+  gPad->SetGridy();
+  
+  ppHist = GetProjections(i, j, 2, &label);
+  
+  for (Int_t centr=0; centr<6; centr++)
+    {
+      if (centr >= 4 || centr == 2)
+       continue;
+      const char* label = 0;
+      TH1* hist = GetProjections(i, j, centr, &label);
+      hist->SetStats(0);
+      hist->SetLineColor(centr+1);
+      hist->Divide(ppHist);
+      hist->GetYaxis()->SetRangeUser(1, 3);
+      hist->Draw((centr == 0) ? "" : "SAME");
+    }
+  legend->Draw();
+  
+  c->SaveAs(Form("phi_%d_%d_ratio.png", i, j));
+}
+
+void CompareProjections(const char* fileName, Int_t i, Int_t j, Int_t centr)
+{
+  TFile::Open(fileName);
+
+  const char* label = 0;
+  TH1* etaHist = 0;
+  TH1* hist = GetProjections(i, j, centr, &label, -1, &etaHist);
+  if (!hist)
+    continue;
+  TH1* hist2 = GetProjections(i, j, centr, &label, 1.2, &etaHist);
+
+  c = new TCanvas;
+  gPad->SetGridx();
+  gPad->SetGridy();
+  
+  hist->SetStats(0);
+  hist->SetMarkerStyle(20);
+  hist->Draw("");
+
+  hist2->SetMarkerStyle(21);
+  hist2->SetMarkerColor(2);
+  hist2->SetLineColor(2);
+  hist2->Draw("SAME");
+}
+
+void DrawProjectionsAll(const char* fileName, Float_t eta = 1, Bool_t etaPhi = kFALSE)
+{
+  DrawProjections(fileName, 0, 1, eta, etaPhi);
+  DrawProjections(fileName, 0, 2, eta, etaPhi);
+  DrawProjections(fileName, 1, 1, eta, etaPhi);
+  DrawProjections(fileName, 1, 2, eta, etaPhi);
+  DrawProjections(fileName, 1, 3, eta, etaPhi);
+  DrawProjections(fileName, 2, 3, eta, etaPhi);
+}
+
+void DrawProjectionsRidge(const char* fileName, Int_t i, Int_t j, Int_t centr1, Int_t centr2)
+{
+  Float_t etaLimit = 1.0;
+  Float_t outerLimit = 1.6;
+  
+  TFile::Open(fileName);
+  
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr1));
+  TH2* hist2 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr2));
+  
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
+
+  SubtractEtaGapNS(hist1, etaLimit, outerLimit, kFALSE);
+  
+  proj1y = ((TH2*) hist1)->ProjectionY("proj1y", hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
+  proj1x = ((TH2*) hist1)->ProjectionX("proj1x", hist1->GetYaxis()->FindBin(-1.79), hist1->GetYaxis()->FindBin(1.79));
+  
+  Float_t etaPhiScale = 1.0 * (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1) / (hist1->GetYaxis()->FindBin(1.79) - hist1->GetYaxis()->FindBin(-1.79) + 1);
+  Printf("%f", etaPhiScale);
+  
+  proj1y->GetXaxis()->SetTitleOffset(1);
+  proj1x->GetXaxis()->SetTitleOffset(1);
+  
+  clone = (TH1*) proj1x->Clone();
+  clone->Fit("pol0", "0", "", 1.2, 1.8);
+  Float_t zyam = clone->GetFunction("pol0")->GetParameter(0);
+
+  proj1x->Add(new TF1("func", "-1", -100, 100), zyam);
+  proj1y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
+  
+  new TCanvas("c", "c", 800, 800);
+  gPad->SetLeftMargin(0.15);
+  //   hist1->SetTitle("");
+
+  hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  hist1->GetXaxis()->SetTitleOffset(1.5);
+  hist1->GetYaxis()->SetTitleOffset(2);
+  hist1->SetStats(kFALSE);
+  hist1->DrawCopy("SURF1");
+  
+  proj2y = ((TH2*) hist2)->ProjectionY("proj2y", hist1->GetXaxis()->FindBin(-0.5), hist1->GetXaxis()->FindBin(0.5));
+  proj2x = ((TH2*) hist2)->ProjectionX("proj2x", hist1->GetYaxis()->FindBin(-1.79), hist1->GetYaxis()->FindBin(1.79));
+
+  //   proj2y->Scale(1.0 / (hist1->GetXaxis()->FindBin(0.5) - hist1->GetXaxis()->FindBin(-0.5) + 1));
+  //   proj2x->Scale(1.0 / (hist1->GetYaxis()->FindBin(1.79) - hist1->GetYaxis()->FindBin(-1.79) + 1));
+  clone = (TH1*) proj2x->Clone();
+  clone->Fit("pol0", "0", "", 1.2, 1.8);
+  zyam = clone->GetFunction("pol0")->GetParameter(0);
+
+  proj2x->Add(new TF1("func", "-1", -100, 100), zyam);
+  proj2y->Add(new TF1("func", "-1", -100, 100), zyam * etaPhiScale);
+
+  proj2y->SetLineColor(2); proj2x->SetLineColor(2);
+  
+  new TCanvas; proj1y->Draw(); proj2y->Draw("SAME"); gPad->SetGridx(); gPad->SetGridy();
+  new TCanvas; proj1x->Draw(); proj2x->Draw("SAME"); gPad->SetGridx(); gPad->SetGridy();
+  
+  new TCanvas("c2", "c2", 800, 800);
+  gPad->SetLeftMargin(0.15);
+  //   hist2->SetTitle("");
+  hist2->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  hist2->GetXaxis()->SetTitleOffset(1.5);
+  hist2->GetYaxis()->SetTitleOffset(2);
+  hist2->SetStats(kFALSE);
+  hist2->DrawCopy("SURF1");
+}
+
+void CalculateIAA(Int_t i, Int_t j, Int_t centr, Float_t etaBegin, Double_t& nsPeak, Double_t& nsPeakE, Double_t& asPeak, Double_t& asPeakE, Double_t& nsRidge, Double_t& nsRidgeE, char** label, Bool_t subtractRidge)
+{
+  TH2* hist = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+  if (!hist)
+    {
+      nsPeak = 0;
+      nsPeakE = 0;
+      asPeak = 0;
+      asPeakE = 0;
+      nsRidge = 0;
+      nsRidgeE = 0;
+      return;
+    }
+  
+  //   new TCanvas; hist->Draw("COLZ");
+  
+  TH2* hist1 = (TH2*) hist->Clone();
+  hist1->Reset();
+  
+  // copy to one quadrant
+  for (Int_t x=1; x<=hist->GetNbinsX(); x++)
+    for (Int_t y=1; y<=hist->GetNbinsY(); y++)
+      {
+       Int_t xTarget = x;
+       if (hist->GetXaxis()->GetBinCenter(x) < 0)
+         xTarget = hist->GetXaxis()->FindBin(-1.0 * hist->GetXaxis()->GetBinCenter(x));
+       else if (hist->GetXaxis()->GetBinCenter(x) > TMath::Pi())
+         xTarget = hist->GetXaxis()->FindBin(TMath::TwoPi() - hist->GetXaxis()->GetBinCenter(x));
+      
+       Int_t yTarget = y;
+       if (hist->GetYaxis()->GetBinCenter(y) < 0)
+         yTarget = hist->GetYaxis()->FindBin(-1.0 * hist->GetYaxis()->GetBinCenter(y));
+      
+       //       Printf("%d %d --> %d %d", x, y, xTarget, yTarget);
+      
+       Float_t value = 0;
+       Float_t error = 0;
+       value += hist->GetBinContent(x, y);
+       error += hist->GetBinError(x, y) * hist->GetBinError(x, y);
+
+       value += hist1->GetBinContent(xTarget, yTarget);
+       error += hist1->GetBinError(xTarget, yTarget) * hist1->GetBinError(xTarget, yTarget);
+
+       error = TMath::Sqrt(error);
+      
+       hist1->SetBinContent(xTarget, yTarget, value);
+       hist1->SetBinError(xTarget, yTarget, error);
+      }
+  
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  new TCanvas; hist1->Draw("COLZ");
+  //   new TCanvas; hist1->Draw("SURF1");
+  
+  tokens = TString(hist1->GetTitle()).Tokenize("-");
+  centralityStr = new TString;
+  if (tokens->GetEntries() > 1)
+    {
+      *centralityStr = tokens->At(0)->GetName();
+      *centralityStr = *centralityStr + "-" + tokens->At(1)->GetName();
+    }
+  *label = centralityStr->Data();
+
+  Int_t phi1 = hist1->GetXaxis()->FindBin(0.0001);
+  Int_t phi3 = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.3);
+  Int_t phi2 = phi3 - 1;
+  Int_t phi4 = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.1);
+  Int_t phi5 = phi4 + 1;
+  Int_t phi6 = hist1->GetXaxis()->FindBin(TMath::Pi() - 0.0001);
+  //   Printf("phi = %d %d %d %d %d %d", phi1, phi2, phi3, phi4, phi5, phi6);
+  
+  Int_t eta1 = hist1->GetYaxis()->FindBin(0.0001);
+  Int_t eta2 = hist1->GetYaxis()->FindBin(etaBegin - 0.0001);
+  Int_t eta3 = eta2 + 1;
+  Int_t eta4 = hist1->GetYaxis()->FindBin(etaMax - 0.0001);
+  //   Printf("eta = %d %d %d %d", eta1, eta2, eta3, eta4);
+  
+  Double_t zyamYieldE, zyamYieldE1, zyamYieldE2;
+  nsPeak              = hist1->IntegralAndError(phi1, phi2, eta1, eta2, nsPeakE, "width");
+  nsRidge             = hist1->IntegralAndError(phi1, phi2, eta3, eta4, nsRidgeE, "width");
+  asPeak              = hist1->IntegralAndError(phi5, phi6, eta1, eta4, asPeakE, "width");
+  Double_t zyamYield  = hist1->IntegralAndError(phi3, phi4, eta1, eta4, zyamYieldE, "width");
+  Double_t zyamYield1 = hist1->IntegralAndError(phi3, phi4, eta1, eta2, zyamYieldE1, "width");
+  Double_t zyamYield2 = hist1->IntegralAndError(phi3, phi4, eta3, eta4, zyamYieldE2, "width");
+  
+  // factor 4 from folding to one quadrant above
+  Double_t zyamYield2Density = zyamYield2 / (hist1->GetXaxis()->GetBinUpEdge(phi4) - hist1->GetXaxis()->GetBinLowEdge(phi3)) / (hist1->GetYaxis()->GetBinUpEdge(eta4) - hist1->GetYaxis()->GetBinLowEdge(eta3)) / 4;
+  Double_t zyamYieldE2Density = zyamYieldE2 / (hist1->GetXaxis()->GetBinUpEdge(phi4) - hist1->GetXaxis()->GetBinLowEdge(phi3)) / (hist1->GetYaxis()->GetBinUpEdge(eta4) - hist1->GetYaxis()->GetBinLowEdge(eta3)) / 4;
+  
+  Double_t nsZyamScaling = 1.0 * (phi2 - phi1 + 1) / (phi4 - phi3 + 1);
+  Double_t asZyamScaling = 1.0 * (phi6 - phi5 + 1) / (phi4 - phi3 + 1);
+  //   Printf("%f %f", nsZyamScaling, asZyamScaling);
+  
+  nsPeak -= zyamYield1 * nsZyamScaling;
+  nsPeakE = TMath::Sqrt(zyamYieldE1 * zyamYieldE1 * nsZyamScaling * nsZyamScaling + nsPeakE * nsPeakE);
+  
+  nsRidge -= zyamYield2 * nsZyamScaling;
+  nsRidgeE = TMath::Sqrt(zyamYieldE2 * zyamYieldE2 * nsZyamScaling * nsZyamScaling + nsRidgeE * nsRidgeE);
+
+  asPeak -= zyamYield * asZyamScaling;
+  asPeakE = TMath::Sqrt(zyamYieldE * zyamYieldE * asZyamScaling * asZyamScaling + asPeakE * asPeakE);
+  
+  if (subtractRidge)
+    {
+      Double_t nsRidgeScaling = 1.0 * (eta2 - eta1 + 1) / (eta4 - eta3 + 1);
+      Double_t asRidgeScaling = 1.0 * (eta4 - eta1 + 1) / (eta4 - eta3 + 1);
+    
+      nsPeak -= nsRidge * nsRidgeScaling;
+      asPeak -= nsRidge * asRidgeScaling;
+      nsPeakE = TMath::Sqrt(nsRidgeE * nsRidgeE * nsRidgeScaling * nsRidgeScaling + nsPeakE * nsPeakE);
+      asPeakE = TMath::Sqrt(nsRidgeE * nsRidgeE * asRidgeScaling * asRidgeScaling + asPeakE * asPeakE);
+    }
+
+  nsRidge /= (etaMax - etaBegin) * 2;
+  nsRidgeE /= (etaMax - etaBegin) * 2;
+  
+  Printf("Peak yields (%d %d %d): %f +- %f; %f +- %f; %f +- %f; %f +- %f", i, j, centr, nsPeak, nsPeakE, asPeak, asPeakE, nsRidge, nsRidgeE, zyamYield2Density, zyamYieldE2Density);
+
+  //   Printf("");
+}
+
+void PlotIAA(const char* fileName)
+{
+  Int_t colors[] = { 1, 2, 3, 4, 6, 7 };
+  Int_t markers[] = { 20, 21, 22, 23, 24, 25 };
+  
+  if (0)
+    {
+      Int_t n = 6;
+      Int_t is[] = { 0, 0, 1, 1, 1, 2 };
+      Int_t js[] = { 1, 2, 1, 2, 3, 3 };
+  
+      Int_t centralityBins = 6;
+      Float_t centralityX[] = { 10, 30, 110, 50, 80, 120 };
+      Float_t centralityEX[] = { 10, 10, 0, 10, 20, 0 };
+    }
+  else if (0)
+    {
+      Int_t n = 3;
+      Int_t is[] = { 0, 1, 2, 3 };
+      Int_t js[] = { 1, 2, 3, 4 };
+
+      Int_t centralityBins = 4;
+      Float_t centralityX[] = { 1.5, 6.5, 30, 75 };
+      Float_t centralityEX[] = { 1.5, 2.5, 20, 25 };
+    }
+  else
+    {
+      Int_t n = 1;
+      Int_t is[] = { 0 };
+      Int_t js[] = { 1 };
+  
+      Int_t centralityBins = 6;
+      Float_t centralityX[] = { 10, 30, 110, 50, 80, 120 };
+      Float_t centralityEX[] = { 10, 10, 0, 10, 20, 0 };
+    }
+  
+  const char* graphTitles[] = { "NS Yield", "AS Yield", "NS Ridge" };
+
+  TCanvas* canvas[3];
+  for (Int_t ci=0; ci<3; ci++)
+    {
+      canvas[ci] = new TCanvas;
+      gPad->SetGridx();
+      gPad->SetGridy();
+    }
+  
+  TFile::Open(fileName);
+
+  legend = new TLegend(0.4, 0.6, 0.99, 0.99);
+  legend->SetFillColor(0);
+
+  for (Int_t i=0; i<n; i++)
+    {
+      TGraphErrors* graph[5];
+      for (Int_t ci=0; ci<5; ci++)
+       graph[ci] = new TGraphErrors;
+
+      char* label = 0;
+      for (Int_t c=0; c<centralityBins; c++)
+       {
+         Double_t nsYield, nsError, asYield, asError, nsRidge, nsRidgeE;
+         CalculateIAA(is[i], js[i], c, 1.2, nsYield, nsError, asYield, asError, nsRidge, nsRidgeE, &label, kTRUE);
+
+         AddPoint(graph[0], centralityX[c], nsYield, centralityEX[c], nsError);
+         AddPoint(graph[1], centralityX[c], asYield, centralityEX[c], asError);
+         AddPoint(graph[2], centralityX[c], nsRidge, centralityEX[c], nsRidgeE);
+
+         //       if (c != 2 && c != 5)
+         {
+           CalculateIAA(is[i], js[i], c, 1.2, nsYield, nsError, asYield, asError, nsRidge, nsRidgeE, &label, kFALSE);
+           AddPoint(graph[3], centralityX[c], nsYield, 0, nsError);
+           AddPoint(graph[4], centralityX[c], asYield, 0, asError);
+         }
+       }
+
+      for (Int_t ci=0; ci<3; ci++)
+       {
+         canvas[ci]->cd();
+         graph[ci]->SetMarkerStyle(markers[i]);
+         graph[ci]->SetMarkerColor(colors[i]);
+         graph[ci]->SetLineColor(colors[i]);
+         graph[ci]->Draw((i == 0) ? "AP" : "PSAME");
+         graph[ci]->GetYaxis()->SetTitle(graphTitles[ci]);
+         graph[ci]->GetYaxis()->SetRangeUser(0, 1);
+       }
+      for (Int_t ci=3; ci<5; ci++)
+       {
+         canvas[ci-3]->cd();
+         graph[ci]->SetLineColor(colors[i]);
+         graph[ci]->Sort();
+         graph[ci]->Draw("LSAME");
+       }
+    
+      legend->AddEntry(graph[0], label, "P");
+    }
+  
+  canvas[0]->cd();
+  legend->Draw();
+  canvas[0]->SaveAs("ns_yield.png");
+
+  canvas[1]->cd();
+  legend->Draw();
+  canvas[1]->SaveAs("as_yield.png");
+
+  canvas[2]->cd();
+  legend->Draw();
+  canvas[2]->SaveAs("ns_ridge.png");
+}
+
+void Draw2D(const char* fileName, Int_t i, Int_t j, Int_t centr)
+{
+  TFile::Open(fileName);
+  
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+  if (!hist1)
+    return 0;
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  
+  //   hist1->Rebin2D(2, 2); hist1->Scale(0.25);
+  
+  hist1->GetYaxis()->SetRangeUser(-1.79, 1.79);
+  new TCanvas; 
+  hist1->DrawCopy("SURF1");
+}
+
+void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr)
+{
+  CreateGraphStructure();
+
+  CorrelationSubtraction(fileName, i, j, centr, graphs[0], 0);
+}
+
+void CorrelationSubtractionAll(const char* fileName, Int_t centr = 0)
+{
+  Int_t n = 6;
+  Int_t is[] = { 0, 1, 1, 2, 2, 2, 3 };
+  Int_t js[] = { 1, 1, 2, 1, 2, 3, 3 };
+
+  CreateGraphStructure();
+
+  for (Int_t i=0; i<n; i++)
+    CorrelationSubtraction(fileName, is[i], js[i], centr, graphs[0], 0);
+}
+
+Int_t gStudySystematic = 0; // 10 = exclusion zone to 0.5; 11 = exclusion off; 12 = exclusion 0.8 and mirror to AS; 13 = scale peripheral; 14 = exclusion 1.2; 20 = non closure; 30 = baseline; 40 = track cuts;  41 = pair cuts; 50/51 = cms comparison; 60 = other peripheral bin ; 70 = 2 Sigma PID
+Int_t gParticleSpecie = 0;//0:Hadrons 1:Pions 2:Kaons 3:Protons
+const char* gParticleSpecieName[] = {"Hadrons","Pions","Kaons","Protons"};
+
+void SetParticleSpecie(Int_t ParticleSpecie=0){
+  gParticleSpecie=ParticleSpecie;
+}
+
+void CorrelationSubtractionHistogram(const char* fileName, Bool_t centralityAxis = kTRUE, const char* outputFileName = "graphs.root",const char* inputFileNameHadrons)
+{
+  if (gBinning == 0)
+  {
+    Int_t n = 6; //merging og 2.5 - 4
+    Int_t is[] =    { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
+    Int_t js[] =    { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
+  }
+  
+  if (gBinning == 1)
+  {
+    // for extraction when trigger pT is from whole range
+    Int_t n = 7;
+    Int_t is[] =    { 0, 0, 0, 0, 0, 0, 0 };
+    Int_t js[] =    { 1, 2, 3, 4, 5, 6, 9 };
+    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1 };
+  }
+
+  if (gBinning == 2)
+  {
+    // default with all bins
+    Int_t n = 15+6+7;
+    // sorted by pT,T
+    Int_t is[] =    { 0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6 };
+    Int_t js[] =    { 1, 1, 2, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 7 };
+    Bool_t symm[] = { 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 };
+  }
+
+  if (gBinning == 3)
+  {
+    // less bins for v3
+    Int_t n = 3;
+    Int_t is[] = { 0, 1, 2 };
+    Int_t js[] = { 1, 2, 3 };
+    Bool_t symm[] = { 1, 1, 1 };
+  }
+
+  if (gBinning == 4)
+  {
+    // default with all but only symmetric bins
+    Int_t n = 6;
+    // sorted by pT,T
+    Int_t is[] =    { 0, 1, 2, 3, 4, 5, 8 };
+    Int_t js[] =    { 1, 2, 3, 4, 5, 6, 9 };
+    Bool_t symm[] = { 1, 1, 1, 1, 1, 1, 1 };
+  }
+
+  Int_t centralityBins =(gDoSubtraction)?4:5;
+
+  Int_t colors[] = { 1, 2, 3, 4, 6, 7 };
+  Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 };
+
+  CreateGraphStructure();
+  
+  Float_t yMax[] = { 0.1, 2, 0.4, 4, 1.5, 50, 50, 4 };
+
+  TCanvas* canvas[8];
+  for (Int_t ci=0; ci<8; ci++)
+    {
+      canvas[ci] = new TCanvas;
+      gPad->SetGridx();
+      gPad->SetGridy();
+      dummy = new TH2F("dummy", Form(";%s;%s", (centralityAxis) ? "Event class" : "<mult>", graphTitles[(ci < 3) ? ci*2 : ci+3]), 100, 0,(gDoSubtraction)? 60:100, 100, 0, yMax[ci]);
+      if (ci == 0)
+       dummy->GetYaxis()->SetTitle("Yield");
+      if (ci == 1)
+       dummy->GetYaxis()->SetTitle("Width");
+      if (ci == 2)
+       dummy->GetYaxis()->SetTitle("v2 , v3");
+      dummy->SetStats(0);
+      dummy->Draw();
+    }
+  
+  legend = new TLegend(0.4, 0.6, 0.99, 0.99);
+  legend->SetFillColor(0);
+
+  legend2 = new TLegend(0.4, 0.6, 0.99, 0.99);
+  legend2->SetFillColor(0);
+
+  for (Int_t i=0; i<n; i++)
+    {
+      for (Int_t c=0; c<centralityBins; c++)
+       {
+         if (c == 2)
+           continue;
+      
+         Double_t nsYield, nsYieldE, asYield, asYieldE, nsWidth, nsWidthE, asWidth, asWidthE, v2, v2E, v3, v3E, centrality;
+         CorrelationSubtraction(fileName, is[i], js[i], c, graphs[i], 1.5 * (-n/2 + i), kTRUE, symm[i], centralityAxis);
+       }
+
+      for (Int_t ci=0; ci<6; ci++)
+       {
+         canvas[ci/2]->cd();
+
+         graphs[i][ci]->SetMarkerStyle(markers[i%14]);
+         graphs[i][ci]->SetMarkerColor((ci % 2 == 0) ? 1 : 2);
+         graphs[i][ci]->SetLineColor((ci % 2 == 0) ? 1 : 2);
+         graphs[i][ci]->GetXaxis()->SetTitle(dummy->GetXaxis()->GetTitle());
+         graphs[i][ci]->Draw("PSAME");
+       }
+    
+      for (Int_t ci=6; ci<11; ci++)
+       {
+         canvas[ci-3]->cd();
+
+         graphs[i][ci]->SetMarkerStyle(markers[i%14]);
+         graphs[i][ci]->SetMarkerColor(1);
+         graphs[i][ci]->SetLineColor(1);
+         graphs[i][ci]->GetXaxis()->SetTitle(dummy->GetXaxis()->GetTitle());
+         graphs[i][ci]->Draw("PSAME");
+       }
+
+      legend->AddEntry(graphs[i][0], graphs[i][0]->GetTitle(), "P");
+      if (symm[i])
+       legend2->AddEntry(graphs[i][0], graphs[i][0]->GetTitle(), "P");
+    }
+  
+  const char* fileNames[] = { "ridge_yield.png", "ridge_width.png", "v2_v3.png", "ridge_yield_ratio.png", "v3_over_v2.png", "remaining_peak.png", "remaining_jet.png", "chi2ndf.png" };
+  
+  for (Int_t ci=0; ci<8; ci++)
+    {
+      canvas[ci]->cd();
+      if (ci == 2 || ci == 4)
+       legend2->Draw();
+      else
+       legend->Draw();
+      canvas[ci]->SaveAs(fileNames[ci]);
+    }
+  
+  if (gStudySystematic != 0)
+    Printf(">>>>>>>>>>>> WARNING: gStudySystematic set to %d", gStudySystematic);
+  if (gParticleSpecie != 0){
+    Printf("\033[1;31m\n>>Unfold v2 of particle %d with hadron file %s\033[m", gParticleSpecie,inputFileNameHadrons);
+    UnfoldUsingVnHadrons(graphs,inputFileNameHadrons,4);  
+    Printf("\033[1;31m\n>>Unfold v3 of particle %d with hadron file %s\033[m", gParticleSpecie,inputFileNameHadrons);
+    UnfoldUsingVnHadrons(graphs,inputFileNameHadrons,5);  
+  }
+  WriteGraphs(outputFileName);
+}
+
+void CorrelationSubtraction(const char* fileName, Int_t i, Int_t j, Int_t centr, TGraphErrors** graph, Float_t axisOffset, Bool_t silent = kFALSE, Bool_t symmetricpT = kFALSE, Bool_t centralityAxis = kTRUE)
+{
+  Printf(Form("Getting dphi_%d_%d_%d from %s", i, j, centr,fileName));
+  TFile::Open(fileName);
+  
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+  TH2* hist2 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, ((gStudySystematic == 60) ? 6 : 4)));
+  TH1* refMult = (TH1*) gFile->Get("refMult");
+
+  if (!hist1 || !hist2)
+    return;
+
+  hist1 = (TH2*) hist1->Clone();
+  hist2 = (TH2*) hist2->Clone();
+  
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+  hist2->Scale(1.0 / hist2->GetYaxis()->GetBinWidth(1));
+  
+  //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 190000); // for 1, 2, 0
+  //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 75000); // for 2, 3, 0
+  //   hist1->Scale(1.0 / 283550); hist2->Scale(1.0 / 100000); // for 2, 2, 0
+  
+  tokens = TString(hist1->GetTitle()).Tokenize("-");
+  if (tokens->GetEntries() > 1)
+    {
+      TString centralityStr;
+      centralityStr = tokens->At(0)->GetName();
+      centralityStr = centralityStr + "-" + tokens->At(1)->GetName();
+      //     Printf("%s", centralityStr.Data());
+      for (Int_t gID=0; gID<NHists; gID++)
+       {
+         TString xTitle = graph[gID]->GetXaxis()->GetTitle();
+         TString yTitle = graph[gID]->GetYaxis()->GetTitle();
+         //       Printf("%s", centralityStr.Data());
+         graph[gID]->SetTitle(centralityStr);
+         graph[gID]->GetXaxis()->SetTitle(xTitle);
+         graph[gID]->GetYaxis()->SetTitle(yTitle);
+         //       Printf("%s %s", yTitle.Data(), graph[gID]->GetYaxis()->GetTitle());
+       }
+      Double_t centrality = -1;
+      if (tokens->GetEntries() == 4)
+       centrality = 0.5 * (TString(tokens->At(3)->GetName()).Atoi() + TString(tokens->At(2)->GetName()).Atoi());
+    }
+  
+//   hist1->GetYaxis()->SetRangeUser(-1.99, 1.99); hist2->GetYaxis()->SetRangeUser(-1.99, 1.99);
+  if (!silent)
+    {
+      new TCanvas; 
+      copy = (TH2*) hist1->DrawCopy("SURF1");
+//       copy->Rebin2D(2, 2);
+//       copy->Scale(0.25);
+//       copy->GetYaxis()->SetRangeUser(-1.59, 1.59);
+    
+      new TCanvas; 
+      copy = (TH2*) hist2->DrawCopy("SURF1");
+//       copy->Rebin2D(2, 2);
+//       copy->Scale(0.25);
+//       copy->GetYaxis()->SetRangeUser(-1.59, 1.59);
+    }
+  
+  if (gStudySystematic == 13)
+    {
+      Float_t factor[10][10][10][4];
+    
+      if (1)
+       {
+         // From Leonardo, Apr 18th 2013
+         factor[0][1][0][0] = 0.130;
+         factor[0][1][1][0] = 0.115;
+         factor[0][1][3][0] = 0.089;
+         factor[1][2][0][0] = 0.106;
+         factor[1][2][1][0] = 0.088;
+         factor[1][2][3][0] = 0.069;
+         factor[2][3][0][0] = 0.085;
+         factor[2][3][1][0] = 0.074;
+         factor[2][3][3][0] = 0.056;
+         factor[3][4][0][0] = 0.055;
+         factor[3][4][1][0] = 0.059;
+         factor[3][4][3][0] = 0.034;
+         factor[4][5][0][0] = 0.021;
+         factor[4][5][1][0] = 0.034;
+         factor[4][5][3][0] = 0.025;
+         factor[0][1][0][1] = 0.134;
+         factor[0][1][1][1] = 0.121;
+         factor[0][1][3][1] = 0.092;
+         factor[1][2][0][1] = 0.110;
+         factor[1][2][1][1] = 0.090;
+         factor[1][2][3][1] = 0.071;
+         factor[2][3][0][1] = 0.085;
+         factor[2][3][1][1] = 0.068;
+         factor[2][3][3][1] = 0.048;
+         factor[3][4][0][1] = 0.055;
+         factor[3][4][1][1] = 0.066;
+         factor[3][4][3][1] = 0.045;
+         factor[4][5][0][1] = -0.010;
+         factor[4][5][1][1] = 0.002;
+         factor[4][5][3][1] = 0.004;
+         factor[0][1][0][2] = 0.040;
+         factor[0][1][1][2] = 0.062;
+         factor[0][1][3][2] = 0.050;
+         factor[1][2][0][2] = 0.108;
+         factor[1][2][1][2] = 0.099;
+         factor[1][2][3][2] = 0.096;
+         factor[2][3][0][2] = 0.083;
+         factor[2][3][1][2] = 0.081;
+         factor[2][3][3][2] = 0.061;
+         factor[3][4][0][2] = 0.049;
+         factor[3][4][1][2] = 0.072;
+         factor[3][4][3][2] = 0.024;
+         factor[4][5][0][2] = 0.023;
+         factor[4][5][1][2] = 0.049;
+         factor[4][5][3][2] = 0.018;
+         factor[0][1][0][3] = -0.192;
+         factor[0][1][1][3] = 0.003;
+         factor[0][1][3][3] = 0.093;
+         factor[1][2][0][3] = 0.001;
+         factor[1][2][1][3] = -0.008;
+         factor[1][2][3][3] = 0.054;
+         factor[2][3][0][3] = 0.083;
+         factor[2][3][1][3] = 0.090;
+         factor[2][3][3][3] = 0.071;
+         factor[3][4][0][3] = 0.140;
+         factor[3][4][1][3] = 0.163;
+         factor[3][4][3][3] = 0.095;
+         factor[4][5][0][3] = 0.078;
+         factor[4][5][1][3] = 0.055;
+         factor[4][5][3][3] = 0.081;
+       }
+      Printf("i=%d  j=%d centr=%d",i,j,centr);
+      Printf("Using factor[i][j][centr][gParticleSpecie] %f", factor[i][j][centr][gParticleSpecie]);
+      hist2->Scale(1.0 + factor[i][j][centr][gParticleSpecie]);
+    }
+  
+  const Float_t etaFlat = 1.2;
+  const Float_t centralEta = 0.5;
+
+  projCentral = hist1->ProjectionY(Form("%s_proj3y", hist1->GetName()), hist1->GetXaxis()->FindBin(-0.7), hist1->GetXaxis()->FindBin(0.7));
+  projCentral->Scale(hist1->GetXaxis()->GetBinWidth(1));
+  projCentral2 = hist2->ProjectionY(Form("%s_proj4y", hist2->GetName()), hist2->GetXaxis()->FindBin(-0.7), hist2->GetXaxis()->FindBin(0.7));
+  projCentral2->Scale(hist2->GetXaxis()->GetBinWidth(1));
+
+  Float_t yEta1 = projCentral->FindBin(-etaMax+0.01);
+  Float_t yEta2 = projCentral->FindBin(-etaFlat-0.01);
+  Float_t yEta3 = projCentral->FindBin(-centralEta+0.01);
+  Float_t yEta4 = projCentral->FindBin(centralEta-0.01);
+  Float_t yEta5 = projCentral->FindBin(etaFlat+0.01);
+  Float_t yEta6 = projCentral->FindBin(etaMax-0.01);
+  Float_t yBaseline = (projCentral->Integral(yEta1, yEta2) + projCentral->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
+  Float_t yPeak = projCentral->Integral(yEta3, yEta4, "width") - yBaseline * (yEta4 - yEta3 + 1) * projCentral->GetBinWidth(1);
+  Printf("Peak (unsubtracted): %f %f", yBaseline, yPeak);
+
+  if (!silent)
+    {
+      Float_t yBaseline = (projCentral->Integral(yEta1, yEta2) + projCentral->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
+      Float_t yBaseline2 = (projCentral2->Integral(yEta1, yEta2) + projCentral2->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
+      Printf("baselines: %f %f", yBaseline, yBaseline2);
+    
+      projCentral2->Add(new TF1("flat", "1", -5, 5), yBaseline - yBaseline2);
+    
+      new TCanvas; 
+      projCentral->DrawCopy();
+      projCentral2->SetLineColor(2);
+      projCentral2->DrawCopy("SAME");
+    }
+  
+  Float_t xValue1 = -1;
+  Float_t xValue2 = -1;
+  Float_t xValue1vn = -1;
+  Float_t xValue2vn = -1;
+  if (centralityAxis)
+    {
+      xValue1 = centrality + axisOffset;
+      xValue2 = xValue1 + 0.5;
+      xValue1vn = centrality + (Int_t) axisOffset;
+      xValue2vn = xValue1vn + 0.5;
+    }
+  else
+    {
+      if (centrality > 0)
+       {
+         xValue1 = refMult->GetBinContent(refMult->GetXaxis()->FindBin(centrality));
+         xValue2 = xValue1;
+         xValue1vn = xValue1;
+         xValue2vn = xValue1;
+         //       Printf("%f %f %f", centrality, xValue1, xValue2);
+       }
+    }
+  
+  Float_t exclusion = (gDoEtaGap)?2.0:0.;
+  if (gStudySystematic == 10)
+    exclusion = 0.5;
+  else if (gStudySystematic == 11)
+    exclusion = 0.0;
+  else if (gStudySystematic == 14 || gStudySystematic == 51)
+    exclusion = 1.2;
+  
+  Int_t eta1 = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(-etaMax + 0.001)));
+  Int_t eta2 = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(-exclusion - 0.001)));
+  Int_t eta3 = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(exclusion + 0.001)));
+  Int_t eta6 = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(etaMax - 0.001)));
+  Int_t phi1Z = hist1->GetXaxis()->FindBin(TMath::Pi()/2 - 0.2);
+  Int_t phi2Z = hist1->GetXaxis()->FindBin(TMath::Pi()/2 + 0.2);
+  
+  Printf("%d %d %d %d", eta1, eta2, eta3, eta6);
+
+  //   new TCanvas; tmpProj = hist1->ProjectionX("tmp", eta1, eta6); tmpProj->Scale(hist1->GetYaxis()->GetBinWidth(1)); tmpProj->Scale(1.0 / (etaMax * 2)); tmpProj->Draw();
+  
+  Double_t baseLineE;
+  Float_t baseLine = hist1->IntegralAndError(phi1Z, phi2Z, eta1, eta6, baseLineE);
+  baseLine /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
+  baseLineE /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
+  Double_t baseLineE2;
+  Float_t baseLine2 = hist2->IntegralAndError(phi1Z, phi2Z, eta1, eta6, baseLineE2);
+  baseLine2 /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
+  baseLineE2 /= (phi2Z - phi1Z + 1) * (eta6 - eta1 + 1);
+  
+  // get baseline of peripheral bin with large eta gap
+  Int_t eta1Alternative = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(-etaMax + 0.001)));
+  Int_t eta2Alternative = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(-1.2 - 0.001)));
+  Int_t eta3Alternative = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(1.2 + 0.001)));
+  Int_t eta6Alternative = TMath::Min(hist1->GetNbinsY(), TMath::Max(1, hist1->GetYaxis()->FindBin(etaMax - 0.001)));
+  Int_t phi1ZAlternative = hist1->GetXaxis()->FindBin(-TMath::Pi()/2 + 0.001);
+  Int_t phi2ZAlternative = hist1->GetXaxis()->FindBin( TMath::Pi()/2 - 0.001);
+  
+  Printf("%d %d %d %d", eta1Alternative, eta2Alternative, eta3Alternative, eta6Alternative);
+  
+  Float_t baseLine2NSGapE = 0, baseLine2NSGapE2 = 0;
+  Float_t baseLine2NSGap = 0;
+  Int_t bins = 0;
+  if (eta1Alternative != eta2Alternative)
+  {
+    baseLine2NSGap += hist2->IntegralAndError(phi1ZAlternative, phi2ZAlternative, eta1Alternative, eta2Alternative, baseLine2NSGapE);
+    bins += eta2Alternative - eta1Alternative + 1;
+  }
+  if (eta3Alternative != eta6Alternative)
+  {
+    baseLine2NSGap += hist2->IntegralAndError(phi1ZAlternative, phi2ZAlternative, eta3Alternative, eta6Alternative, baseLine2NSGapE2);
+    bins += eta6Alternative - eta3Alternative + 1;
+  }
+  baseLine2NSGapE = TMath::Sqrt(baseLine2NSGapE * baseLine2NSGapE + baseLine2NSGapE2 * baseLine2NSGapE2);
+  baseLine2NSGap /= (phi2ZAlternative - phi1ZAlternative + 1) * bins;
+  baseLine2NSGapE /= (phi2ZAlternative - phi1ZAlternative + 1) * bins;
+  
+  if(gDoSubtraction)hist1->Add(hist2, -1);
+  
+  // phi projection
+  // NS
+  proj = hist1->ProjectionX(Form("%s_proj1x", hist1->GetName()));
+
+  //   proj2 = hist1->ProjectionX(Form("%s_proj2x", hist1->GetName()), eta3, eta6);
+//   proj->Add(proj2, 1);
+  
+  // AS
+  if (gDoEtaGapAS)
+  {
+    projAS = hist1->ProjectionX(Form("%s_proj3x", hist1->GetName()), eta1, eta2);
+    projAS2 = hist1->ProjectionX(Form("%s_proj4x", hist1->GetName()), eta3, eta6);
+    projAS->Add(projAS2, 1);
+    
+    projAS->Scale(1.0 * (eta6 - eta1 + 1) / (eta6 - eta3 + 1 + eta2 - eta1 + 1));
+  }
+  else
+  {
+    projAS = hist1->ProjectionX(Form("%s_proj3x", hist1->GetName()), eta1, eta6);
+  }
+  
+  // match NS and AS yield
+//   proj->Scale(1.0 * (eta6 - eta1 + 1) / (eta6 - eta3 + 1 + eta2 - eta1 + 1));
+
+  // copy AS
+//   for (Int_t bin=proj->FindBin(TMath::Pi()/2+0.001); bin<=proj->GetNbinsX(); bin++)
+//     {
+//       //     Printf("%d %f", bin, projAS->GetBinContent(bin));
+//       proj->SetBinContent(bin, projAS->GetBinContent(bin));
+//       proj->SetBinError(bin, projAS->GetBinError(bin));
+//     }
+  
+  if (gStudySystematic == 12)
+    {
+      // get difference between exclusion 0.8 and 0.0, shift to AS, remove from AS
+      Printf("\nParticleSpecie: %d",gParticleSpecie);
+      projAll = hist1->ProjectionX(Form("%s_proj4x", hist1->GetName()), eta2+1, eta3-1);
+      projAll->Add(proj, -1.0 * (eta3-1 - (eta2+1) + 1) / (eta6 - eta1 + 1));
+  
+      if (!silent)
+       {
+         new TCanvas; projAll->DrawCopy();
+       }
+
+      // From Leonardo, Apr 18th 2013
+      Float_t ratioNSAS[10][10][10][4];
+      ratioNSAS[0][1][4][0] = 3.166786;
+      ratioNSAS[1][2][4][0] = 1.595074;
+      ratioNSAS[2][3][4][0] = 1.451555;
+      ratioNSAS[3][4][4][0] = 1.371433;
+      ratioNSAS[4][5][4][0] = 1.581739;
+      ratioNSAS[0][1][4][1] = 3.177672;
+      ratioNSAS[1][2][4][1] = 1.628396;
+      ratioNSAS[2][3][4][1] = 1.497103;
+      ratioNSAS[3][4][4][1] = 1.430801;
+      ratioNSAS[4][5][4][1] = 1.627891;
+      ratioNSAS[0][1][4][2] = 2.820564;
+      ratioNSAS[1][2][4][2] = 1.526886;
+      ratioNSAS[2][3][4][2] = 1.415753;
+      ratioNSAS[3][4][4][2] = 1.195244;
+      ratioNSAS[4][5][4][2] = 1.519867;
+      ratioNSAS[0][1][4][3] = 1.825534;
+      ratioNSAS[1][2][4][3] = 1.143225;
+      ratioNSAS[2][3][4][3] = 1.090177;
+      ratioNSAS[3][4][4][3] = 1.239809;
+      ratioNSAS[4][5][4][3] = 1.091130;
+      
+      Printf("i=%d  j=%d centr=%d",i,j,centr);
+      Printf("Using NS/AS ratio %f", ratioNSAS[i][j][4][gParticleSpecie]);
+
+      for (Int_t k=0; k<=projAll->GetNbinsX()/2; k++)
+       {
+         //     Printf("%d -> %d", i, projAll->GetNbinsX()+1-i);
+         projAll->SetBinContent(projAll->GetNbinsX()+1-k, projAll->GetBinContent(k) / ratioNSAS[i][j][4][gParticleSpecie]);
+         projAll->SetBinError(projAll->GetNbinsX()+1-k, projAll->GetBinError(k) / ratioNSAS[i][j][4][gParticleSpecie]);
+      
+         projAll->SetBinContent(k, 0);
+         projAll->SetBinError(k, 0);
+       }
+      if (!silent)
+       {
+         new TCanvas; projAll->DrawCopy();
+       }
+    
+      //     new TCanvas; proj->DrawCopy();
+      proj->Add(projAll, -1);
+      //     new TCanvas; proj->DrawCopy();
+    }
+  
+//   proj->Scale(hist1->GetYaxis()->GetBinWidth(1));
+//   proj->Scale(1.0 / (etaMax * 2));
+  proj->Scale(1.0 / hist1->GetNbinsY());
+  proj->Rebin(2); proj->Scale(0.5);
+  
+  if(gExcludeNearSidePeakFromPhiProj)
+    {
+      Double_t phicut=1.;
+      for(Int_t ibin=1;ibin<=proj->GetNbinsX();ibin++)
+       if(TMath::Abs(proj->GetBinCenter(ibin))<phicut)
+         {
+           proj->SetBinContent(ibin,0);
+           proj->SetBinError(ibin,0);
+         }
+    }
+  
+  if (gStudySystematic == 20)
+    {
+      Printf(">>>>>>>>>>>> Applying non-closure systematics <<<<<<<<<<<<");
+      Printf(">> Non closure ratio from file non_closure_%s.root",gParticleSpecieName[gParticleSpecie]);
+      file2 = TFile::Open(Form("non_closure_files/non_closure_%s.root",gParticleSpecieName[gParticleSpecie]));
+      non_closure = (TH1*) file2->Get(Form("non_closure_%d_%d_%d", i, j, 0));
+      for (Int_t bin=1; bin<=non_closure->GetNbinsX(); bin++)
+       non_closure->SetBinError(bin, 0);
+    
+      proj->Multiply(non_closure);
+    }
+  
+  //   proj = hist1->ProjectionX("proj1x", hist1->GetYaxis()->FindBin(0.5), eta6);
+
+  // eta projection
+  Int_t binEtaProj1 = 1;
+  Int_t binEtaProj2 = hist1->GetXaxis()->FindBin(-TMath::Pi()/3-0.001);
+  Int_t binEtaProj3 = binEtaProj2 + 1;
+  Int_t binEtaProj4 = hist1->GetXaxis()->FindBin(TMath::Pi()/3-0.001);
+  Int_t binEtaProj5 = binEtaProj4+1;
+  Int_t binEtaProj6 = hist1->GetXaxis()->FindBin(TMath::Pi()-TMath::Pi()/3-0.001);
+  Int_t binEtaProj7 = binEtaProj6+1;
+  Int_t binEtaProj8 = hist1->GetXaxis()->FindBin(TMath::Pi()+TMath::Pi()/3-0.001);
+  Int_t binEtaProj9 = binEtaProj8+1;
+  Int_t binEtaProj10 = hist1->GetNbinsX();
+  Printf("%d %d %d %d %d %d %d %d %d %d", binEtaProj1, binEtaProj2, binEtaProj3, binEtaProj4, binEtaProj5, binEtaProj6, binEtaProj7, binEtaProj8, binEtaProj9, binEtaProj10);
+  
+  etaProj = hist1->ProjectionY(Form("%s_proj1y", hist1->GetName()), binEtaProj3, binEtaProj4);
+  etaProj->Scale(1.0 / (binEtaProj4 - binEtaProj3 + 1));
+  etaProj2 = hist1->ProjectionY(Form("%s_proj2y", hist1->GetName()), binEtaProj7, binEtaProj8);
+  etaProj2->Scale(1.0 / (binEtaProj8 - binEtaProj7 + 1));
+  etaProj3 = hist1->ProjectionY(Form("%s_proj5y", hist1->GetName()), binEtaProj1, binEtaProj2);
+  etaProj3b = hist1->ProjectionY(Form("%s_proj5by", hist1->GetName()), binEtaProj5, binEtaProj6);
+  etaProj3c = hist1->ProjectionY(Form("%s_proj5cy", hist1->GetName()), binEtaProj9, binEtaProj10);
+  etaProj3->Add(etaProj3b);
+  etaProj3->Add(etaProj3c);
+  etaProj3->Scale(1.0 / (binEtaProj2 - binEtaProj1 + 1 + binEtaProj6 - binEtaProj5 + 1 + binEtaProj10 - binEtaProj9 + 1));
+
+  Float_t ySubBaseline = (etaProj->Integral(yEta1, yEta2) + etaProj->Integral(yEta5, yEta6)) / (yEta2 - yEta1 + 1 + yEta6 - yEta5 + 1);
+  Float_t ySubPeak = etaProj->Integral(yEta3, yEta4, "width") - ySubBaseline * (yEta4 - yEta3 + 1) * etaProj->GetBinWidth(1);
+  Printf("Peak (subtracted): %f %f (%.2f%% of unsubtracted peak)", ySubBaseline, ySubPeak, ySubPeak / yPeak * 100);
+  Printf("    factor[%d][%d][%d][%d] = %.3f;", i, j, centr, gParticleSpecie, ySubPeak / yPeak);
+  AddPoint(graph[8], xValue1, 100.0 * ySubPeak / yPeak, 0, 0);
+  
+//   hist1->Rebin2D(4, 4); hist1->Scale(1.0 / 16);
+//   hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
+  
+  TString fitOption = "0I";
+  if (!silent)
+    {
+      c = new TCanvas("c", "c", 600, 600);
+
+      PadFor2DCorr();
+    
+//       hist1->GetYaxis()->SetRangeUser(-1.59, 1.59);
+      hist1->GetXaxis()->SetTitleOffset(1.7);
+      hist1->GetYaxis()->SetTitleOffset(2);
+      hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
+      hist1->GetZaxis()->SetTitleOffset(2);
+      hist1->GetZaxis()->SetNdivisions(504);
+      hist1->SetStats(kFALSE);
+      hist1->GetZaxis()->CenterTitle(kTRUE);
+      hist1->GetYaxis()->CenterTitle(kTRUE);
+      hist1->GetXaxis()->CenterTitle(kTRUE);
+      hist1->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+      hist1->GetYaxis()->SetNdivisions(505);
+    
+      TString label(hist1->GetTitle());
+      hist1->SetTitle("");
+      //label.ReplaceAll("00", "");
+      //label.ReplaceAll(".0", "");
+      TObjArray* objArray = label.Tokenize("-");
+      TPaveText* paveText = new TPaveText(0.03, 0.86, 0.44, 0.97, "BRNDC");
+      paveText->SetTextSize(0.035);
+      paveText->SetFillColor(0);
+      paveText->SetShadowColor(0);
+      paveText->SetBorderSize(0);
+      paveText->SetFillStyle(0);
+      TString tmpStr(objArray->At(0)->GetName());
+      tmpStr.ReplaceAll("p_", "#it{p}_");
+      paveText->AddText(tmpStr + "GeV/#it{c}");
+
+      TString tmpStr(objArray->At(1)->GetName());
+      tmpStr.ReplaceAll("p_", "#it{p}_");
+      paveText->AddText(tmpStr + "GeV/#it{c}");
+  
+      TString label2(hist2->GetTitle());
+      TObjArray* objArray2 = label2.Tokenize("-");
+
+      TPaveText* paveText2 = new TPaveText(0.65, 0.86, 0.98, 0.97, "BRNDC");
+      paveText2->SetTextSize(0.035);
+      paveText2->SetBorderSize(0);
+      paveText2->SetFillColor(0);
+      paveText2->SetFillStyle(0);
+      paveText2->SetShadowColor(0);
+      paveText2->AddText("p-Pb #sqrt{s_{NN}} = 5.02 TeV");
+      if (objArray->GetEntries() > 3 && objArray2->GetEntries() > 3)
+       {
+         TString centrBeginStr;
+         centrBeginStr = objArray->At(2)->GetName();
+         centrBeginStr.ReplaceAll(" ", "");
+         TString centrEndStr;
+         centrEndStr = objArray2->At(2)->GetName();
+         centrEndStr.ReplaceAll(" ", "");
+         if(gDoSubtraction)paveText2->AddText(Form("(%s-%s) - (%s-%s)", centrBeginStr.Data(), objArray->At(3)->GetName(), centrEndStr.Data(), objArray2->At(3)->GetName()));
+         else paveText2->AddText(Form("(%s-%s)", centrBeginStr.Data(), objArray->At(3)->GetName()));
+       }
+      
+      hist1->DrawCopy("SURF1");
+      paveText->Draw();
+      paveText2->Draw();
+      gPad->GetCanvas()->SaveAs(Form("ridge_%d_%d.png", i, j));
+      //gPad->GetCanvas()->SaveAs(Form("ridge_%d_%d.eps", i, j));
+      //gPad->GetCanvas()->SaveAs("fig3a.eps");
+    
+      Float_t fontSize = 0.05;
+
+      c3 = new TCanvas("c3", "c3", 600, 400);
+      gPad->SetLeftMargin(0.12);
+      gPad->SetRightMargin(0.01);
+      gPad->SetBottomMargin(0.12);
+      gPad->SetTopMargin(0.01);
+      etaProj->SetTitle();
+      //     etaProj->GetYaxis()->SetNdivisions(505);
+      etaProj->GetYaxis()->SetLabelSize(fontSize);
+      etaProj->GetXaxis()->SetLabelSize(fontSize);
+      etaProj->GetXaxis()->SetTitleSize(fontSize);
+      etaProj->GetYaxis()->SetTitleSize(fontSize);
+      etaProj->GetYaxis()->SetTitle(kProjYieldTitleEta);
+      etaProj->GetYaxis()->SetTitleOffset(1.1);
+      etaProj->SetStats(0);
+      etaProj->GetXaxis()->SetNdivisions(505);
+      etaProj->Rebin(2); etaProj->Scale(0.5);
+      etaProj2->Rebin(2); etaProj2->Scale(0.5);
+      etaProj3->Rebin(2); etaProj3->Scale(0.5);
+      etaProj->Rebin(2); etaProj->Scale(0.5);
+      etaProj2->Rebin(2); etaProj2->Scale(0.5);
+      etaProj3->Rebin(2); etaProj3->Scale(0.5);
+//       etaProj->GetXaxis()->SetRangeUser(-1.99, 1.99);
+      //     etaProj->GetYaxis()->SetRangeUser(TMath::Min(etaProj->GetMinimum(), etaProj3->GetMinimum()) * 0.8, TMath::Max(etaProj->GetMaximum(), etaProj3->GetMaximum()) * 1.2);
+      etaProj->GetYaxis()->SetRangeUser(proj->GetMinimum(0.00001) * 0.98, proj->GetMaximum() * 1.1);
+      etaProj2->SetLineColor(2); etaProj2->SetMarkerColor(2);
+      etaProj3->SetLineColor(4); etaProj3->SetMarkerColor(4);
+      etaProj->SetMarkerStyle(24);
+      etaProj2->SetMarkerStyle(25);
+      etaProj3->SetMarkerStyle(26);
+      etaProj->DrawCopy("E0 X0");
+      etaProj2->DrawCopy("E0 X0 SAME");
+      etaProj3->DrawCopy("E0 X0 SAME");
+    
+      legend5 = new TLegend(0.48, 0.74, 0.92, 0.97);
+      legend5->SetBorderSize(0);
+      legend5->SetTextSize(fontSize);
+      legend5->SetFillColor(0);
+      legend5->AddEntry(etaProj,  "|#Delta#varphi| < #pi/3", "P");
+      legend5->AddEntry(etaProj2, "|#Delta#varphi - #pi| < #pi/3", "P");
+      //     legend5->AddEntry(etaProj3, "#pi/3 < |#Delta#varphi| < 2#pi/3, #Delta#varphi > 4#pi/3", "P");
+      legend5->AddEntry(etaProj3, "Remaining #Delta#varphi", "P");
+      legend5->Draw();
+    
+      c = new TCanvas("c2", "c2", 600, 400);
+      gPad->SetLeftMargin(0.12);
+      gPad->SetRightMargin(0.01);
+      gPad->SetBottomMargin(0.12);
+      gPad->SetTopMargin(0.01);
+    
+      proj->SetStats(0);
+      //     proj->GetYaxis()->SetNdivisions(505);
+      proj->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+      proj->GetYaxis()->SetTitle(kProjYieldTitlePhi);
+      proj->GetYaxis()->SetTitleOffset(1.1);
+      proj->GetYaxis()->SetLabelSize(fontSize);
+      proj->GetXaxis()->SetLabelSize(fontSize);
+      proj->GetXaxis()->SetTitleSize(fontSize);
+      proj->GetYaxis()->SetTitleSize(fontSize);
+      proj->SetTitle();
+      proj->SetMarkerStyle(21);
+      proj->SetMarkerSize(0.7);
+      //proj->SetMinimum(proj->GetMinimum() * 0.98); proj->SetMaximum(proj->GetMaximum() * 1.065);
+      //proj->GetYaxis()->SetRangeUser(proj->GetMinimum(0.00001) * 0.98,proj->GetMaximum() * 1.6);
+      proj->GetYaxis()->SetRangeUser(proj->GetMinimum(0.00001) * 0.98,proj->GetMaximum() * 1.05);
+      proj->Draw("E0 X0");
+      fitOption = "I";
+    
+      if (0)
+       {
+         Printf("\nPer-trigger yield per unit of delta eta (y) as function of delta phi (x) in the bin %s", centralityStr.Data());
+         Printf("Systematic uncertainties are mostly correlated and affect the baseline. Uncorrelated uncertainties are less than 1%% and not indicated explicitly in the table below.");
+
+         for (Int_t k=1; k<=proj->GetNbinsX(); k++)
+           Printf("x = %.2f, y = %.4f +- %.4f (stat)", proj->GetXaxis()->GetBinCenter(k), proj->GetBinContent(k), proj->GetBinError(k));
+       }
+    }
+  
+  fileProj = TFile::Open("dphi_proj.root", "UPDATE");
+  proj->Write(Form("proj_%d_%d_%d_%d", i, j, centr, gStudySystematic));
+  fileProj->Close();
+
+  TF1* v2 = new TF1("func", "[0]+2*[1]*cos(2*x)", -5, 5);
+  v2->SetLineStyle(2);
+  //   v2->SetLineWidth(1);
+  
+  TF1* v2v3 = new TF1("func", "[0]+2*[1]*cos(2*x)+2*[2]*cos(3*x)+2*[3]*cos(x)", -5, 5);
+
+  v2v3->SetLineColor(1);
+  //   v2v3->FixParameter(2, 0);
+//   proj->Fit(v2, fitOption);
+  //   return;
+  
+  fitOption += "+";
+  proj->Fit(v2v3, fitOption, "E0 X0");
+  
+  if(gDoSubtraction){
+    Float_t min = v2v3->GetMinimum();
+    Printf("Minimum for the yield extraction from the fit %f",min);
+  }else{
+    //get the minimum value (not using GetMinimum() because it can be that same bins are 0, when we remove the NS peak in the delta_phi projection)
+    Float_t min = 99999;
+    for(Int_t ibin=1;ibin<=proj->GetNbinsX();ibin++){
+      if(proj->GetBinContent(ibin)==0)continue;
+      if(proj->GetBinContent(ibin)<min)min=proj->GetBinContent(ibin);
+    }
+    Printf("Minimum for the yield extraction from the projection %f",min);
+  }
+  
+  Float_t diffMinParam0 = v2v3->GetParameter(0) - v2v3->GetMinimum();
+  Printf("Chi2: %f ndf: %d chi2/ndf %f:Min: %f %f", v2v3->GetChisquare(), v2v3->GetNDF(), v2v3->GetChisquare() / v2v3->GetNDF(), min, diffMinParam0);
+  AddPoint(graph[10], xValue1, v2v3->GetChisquare() / v2v3->GetNDF(), 0, 0);
+  
+  AddPoint(graph[11], xValue1, baseLine + diffMinParam0, 0, baseLineE);
+  //   AddPoint(graph[11], xValue1, baseLine, 0, baseLineE);
+
+  
+  
+  Float_t v2value = 0, v2E = 0, v3 = 0, v3E = 0;
+  if(gUseAnalyticalVn){//analytic calculation of vv
+    Printf("vn calculated analytically!");
+    v2value = CalculateVnAnalytically(proj,2,(gDoSubtraction)?baseLine2:0.);//baseline2 is the peripheral one
+    v2E = CalculateVnErrAnalytically(proj,2,(gDoSubtraction)?baseLine2:0.,(gDoSubtraction)?baseLineE2:0.);
+    if (symmetricpT)
+      AddPoint(graph[4], xValue1vn, v2value, 0, v2E);
+    v3 = CalculateVnAnalytically(proj,3,(gDoSubtraction)?baseLine2:0.);//baseline2 is the peripheral one
+    v3E = CalculateVnErrAnalytically(proj,3,(gDoSubtraction)?baseLine2:0.,(gDoSubtraction)?baseLineE2:0);
+    if (symmetricpT)
+      AddPoint(graph[5], xValue2vn, v3, 0, v3E);
+    
+    if (0 && !silent)
+    {
+      proj->Fit("pol0", "+W");
+      Float_t v0 = proj->GetFunction("pol0")->GetParameter(0);
+      Float_t v1 = CalculateVnAnalyticallyRaw(proj,1,(gDoSubtraction)?baseLine2:0.);//baseline2 is the peripheral one
+      Printf("%f", v1);
+      Float_t v4 = 0;//CalculateVnAnalyticallyRaw(proj,4,(gDoSubtraction)?baseLine2:0.);//baseline2 is the peripheral one
+      
+      Float_t scaleBaseline = 1;
+      if (gDoSubtraction)
+       scaleBaseline = v0 * (v0+baseLine2);
+
+      TF1* v1v2v3v4 = new TF1("func", "[0]*(1+2*[1]*cos(2*x)+2*[2]*cos(3*x)+2*[3]*cos(x)+2*[4]*cos(4*x))", -5, 5);
+      v1v2v3v4->SetParameters((Double_t) v0, (Double_t) v2value*v2value / scaleBaseline, (Double_t) v3*v3 * scaleBaseline, (Double_t) v1 / scaleBaseline, (Double_t) v4 / scaleBaseline);
+      v1v2v3v4->SetLineColor(4);
+      v1v2v3v4->Draw("SAME");
+
+      TF1* v2only = (TF1*) v1v2v3v4->Clone("v2only");
+      v2only->SetParameters((Double_t) v0, (Double_t) v2value*v2value);
+      v2only->SetLineColor(4);
+      v2only->SetLineStyle(2);      
+      v2only->Draw("SAME");
+    }
+  }else{//v2 calculated from the fit 
+    Printf("vn calculated from fit!");
+    if (v2v3->GetParameter(1) > 0)
+      {
+//     if (0)
+       {
+         v2value = TMath::Sqrt(v2v3->GetParameter(1) / (baseLine + diffMinParam0));
+         v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + baseLineE * baseLineE / baseLine / baseLine);
+         Printf("-> v2 = %f +- %f (%f %f = %f)", v2value, v2E, baseLine, diffMinParam0, baseLine + diffMinParam0);
+       }
+//     else
+       {
+         v2value = TMath::Sqrt(v2v3->GetParameter(1) / (v2v3->GetParameter(0) + baseLine2NSGap));
+         v2E = 0.5 * v2value * TMath::Sqrt(v2v3->GetParError(1) * v2v3->GetParError(1) / v2v3->GetParameter(1) / v2v3->GetParameter(1) + v2v3->GetParError(0) * v2v3->GetParError(0) / v2v3->GetParameter(0) / v2v3->GetParameter(0) + baseLine2NSGapE * baseLine2NSGapE / baseLine2NSGap / baseLine2NSGap);
+         Printf("-> v2 = %f +- %f (%f %f = %f)", v2value, v2E, v2v3->GetParameter(0), baseLine2NSGap, v2v3->GetParameter(0) + baseLine2NSGap);
+       }
+       if (symmetricpT)
+         AddPoint(graph[4], xValue1vn, v2value, 0, v2E);
+      }
+    
+    if (v2v3->GetParameter(2) > 0)
+      {
+       v3 = TMath::Sqrt(v2v3->GetParameter(2) / (baseLine + diffMinParam0));
+       v3E = 0.5 * v3 * TMath::Sqrt(v2v3->GetParError(2) * v2v3->GetParError(2) / v2v3->GetParameter(2) / v2v3->GetParameter(2) + baseLineE * baseLineE / baseLine / baseLine);
+       if (symmetricpT)
+         AddPoint(graph[5], xValue2vn, v3, 0, v3E);
+      }
+      
+//     Float_t v1 = TMath::Sqrt(TMath::Abs(v2v3->GetParameter(3)) / (baseLine + diffMinParam0));
+//     Float_t v1E = 0.5 * v1 * TMath::Sqrt(v2v3->GetParError(3) * v2v3->GetParError(3) / v2v3->GetParameter(3) / v2v3->GetParameter(3) + baseLineE * baseLineE / baseLine / baseLine);
+//     if (v2v3->GetParameter(3) < 0)
+//       v1 = -v1;
+
+       Float_t v1 = TMath::Abs(v2v3->GetParameter(3)) / (baseLine + diffMinParam0);
+       Float_t v1E = v1 * TMath::Sqrt(v2v3->GetParError(3) * v2v3->GetParError(3) / v2v3->GetParameter(3) / v2v3->GetParameter(3) + baseLineE * baseLineE / baseLine / baseLine);
+
+       if (symmetricpT)
+         AddPoint(graph[12], xValue2vn, v1, 0, v1E);
+//     Printf("v1: %f %f", v1, v1E);
+  }
+  
+  if (v2value > 0 && v3 > 0 && symmetricpT)
+    {
+      Float_t v3v2Ratio = v3 / v2value;
+      Float_t v3v2RatioE = v3v2Ratio * TMath::Sqrt( v2E * v2E / v2value / v2value + v3E * v3E / v3 / v3 );
+      AddPoint(graph[7], xValue1vn, v3v2Ratio, 0, v3v2RatioE);
+    }
+  
+  Printf("Baseline: %f +- %f; v2 = %f +- %f; v3 = %f +- %f", baseLine, baseLineE, v2value, v2E, v3, v3E);
+  
+  if (gStudySystematic == 30)
+    {
+      // alternative way for baseline
+    
+      parabola = new TF1("parabola", "[0] + [1]*(x - [2])**2",  0, TMath::Pi());
+      parabola->SetLineColor(3);
+      parabola->SetParameters(min, -0.1, TMath::Pi() / 2);
+      //     parabola->SetParLimits(1, 0, 1);
+      proj->Fit(parabola, fitOption, "", TMath::Pi() / 2 - 1, TMath::Pi() / 2 + 1);
+      //     proj->Fit(parabola, fitOption, "", 0.1, 2);
+    
+      min = parabola->GetParameter(0);
+      Printf("New baseline: %f", min);
+    }
+
+  fitOption += "R";
+  
+  if (0)
+    {
+      // single gauss fit
+    
+      TF1* gaus1 = new TF1("gaus1", "[0]+gaus(1)", -0.5 * TMath::Pi(), 0.5 * TMath::Pi());
+      gaus1->SetParameters(min, 1, 0, 1);
+      //   gaus1->FixParameter(0, min);
+      gaus1->FixParameter(2, 0);
+      gaus1->SetParLimits(1, 0.001, 10);
+      gaus1->SetParLimits(3, 0.1, 2);
+      gaus1->SetLineColor(3);
+      proj->Fit(gaus1, fitOption);
+    
+      TF1* gaus2 = new TF1("gaus2", "[0]+gaus(1)", 0.5 * TMath::Pi(), 1.5 * TMath::Pi());
+      gaus2->SetParameters(min, 1, TMath::Pi(), 1);
+      gaus2->FixParameter(2, TMath::Pi());
+      gaus2->SetParLimits(1, 0.001, 10);
+      gaus2->SetParLimits(3, 0.1, 2);
+      gaus2->SetLineColor(3);
+      proj->Fit(gaus2, fitOption);
+  
+      AddPoint(graph[2], xValue1, gaus1->GetParameter(3), 0, gaus1->GetParError(3));
+      AddPoint(graph[3], xValue2, gaus2->GetParameter(3), 0, gaus2->GetParError(3));
+    }
+  else if (0)
+    {
+      // combined gauss fit
+    
+      TF1* gausBoth = new TF1("gaus2", "[0]+[1]*(exp(-0.5*(x/[2])**2)+exp(-0.5*((x-TMath::TwoPi())/[2])**2))+[3]*(exp(-0.5*((x-TMath::Pi())/[4])**2)+exp(-0.5*((x+TMath::Pi())/[4])**2))", -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
+      gausBoth->SetParameters(min, 1, 1, 1, 1);
+      //   gausBoth->FixParameter(0, min);
+      gausBoth->SetParLimits(1, 0.0001, 10);
+      gausBoth->SetParLimits(3, 0.0001, 10);
+      gausBoth->SetParLimits(2, 0.1, 4);
+      gausBoth->SetParLimits(4, 0.1, 4);
+      gausBoth->SetLineColor(6);
+    
+      proj->Fit(gausBoth, fitOption);
+  
+      AddPoint(graph[2], xValue1, gausBoth->GetParameter(2), 0, gausBoth->GetParError(3));
+      AddPoint(graph[3], xValue2, gausBoth->GetParameter(4), 0, gausBoth->GetParError(4));
+    }
+  else
+    {
+      // RMS without baseline
+      projRMS = (TH1*) proj->Clone();
+      projRMS->Add(new TF1("f", "1", -10, 10), -min);
+      projRMS->GetXaxis()->SetRangeUser(-TMath::Pi()/2 + 0.01, TMath::Pi()/2 - 0.01);
+      AddPoint(graph[2], xValue1, projRMS->GetRMS(), 0, projRMS->GetRMSError());
+      projRMS->GetXaxis()->SetRangeUser(TMath::Pi()/2 + 0.01, 3 * TMath::Pi()/2 - 0.01);
+      AddPoint(graph[3], xValue2, projRMS->GetRMS(), 0, projRMS->GetRMSError());
+      //     new TCanvas; projRMS->Draw(); return;
+    }
+  
+  //   new TCanvas; gausBoth->Draw(); return;
+
+  if (!silent)
+    {
+      TF1* v2v3_v1 = new TF1("func", "[0]+2*[1]*cos(x)", -5, 5);
+      v2v3_v1->SetParameters(v2v3->GetParameter(0), v2v3->GetParameter(3));
+      v2v3_v1->SetLineStyle(3);
+      v2v3_v1->SetLineColor(4);
+//       v2v3_v1->SetLineWidth(1);
+      v2v3_v1->Draw("SAME");
+      TF1* v2v3_v2 = new TF1("func", "[0]+2*[1]*cos(2*x)", -5, 5);
+      v2v3_v2->SetParameters(v2v3->GetParameter(0), v2v3->GetParameter(1));
+      v2v3_v2->SetLineStyle(2);
+      v2v3_v2->SetLineColor(2);
+//       v2v3_v2->SetLineWidth(1);
+      v2v3_v2->Draw("SAME");
+      TF1* v2v3_v3 = new TF1("func", "[0]+2*[1]*cos(3*x)", -5, 5);
+      v2v3_v3->SetParameters(v2v3->GetParameter(0), v2v3->GetParameter(2));
+      v2v3_v3->SetLineStyle(4);
+      v2v3_v3->SetLineColor(6);
+//       v2v3_v3->SetLineWidth(1);
+      v2v3_v3->Draw("SAME");
+    
+      line = new TLine(-0.5 * TMath::Pi(), min, 1.5 * TMath::Pi(), min);
+      line->SetLineWidth(2);
+      line->SetLineColor(4);
+      line->Draw();
+
+      legend = new TLegend(0.47, 0.65, 0.88, 0.95);
+      legend->SetFillColor(0);
+      legend->SetBorderSize(0);
+
+      legend->AddEntry(proj, "Data", "P");
+      legend->AddEntry(v2v3, "a_{0} + a_{2} cos(2#Delta#varphi) + a_{3} cos(3#Delta#varphi)", "L");
+      legend->AddEntry(v2, "a_{0} + a_{2} cos(2#Delta#varphi)", "L");
+      legend->AddEntry(line, "Baseline for yield extraction", "L");
+
+      legend->SetTextSize(fontSize);
+      legend->Draw();
+    
+      //     paveText3 = (TPaveText*) paveText2->Clone();
+      //     paveText3->SetX1NDC(0.16);
+      //     paveText3->SetY1NDC(0.77);
+      //     paveText3->SetX2NDC(0.42);
+      //     paveText3->SetY2NDC(0.94);
+      //     paveText3->Draw();
+      paveText4 = (TPaveText*) paveText2->Clone();
+      paveText4->SetTextSize(fontSize);
+      paveText4->SetX1NDC(0.16);
+      paveText4->SetY1NDC(0.68);
+      paveText4->SetX2NDC(0.42);
+      paveText4->SetY2NDC(0.96);
+
+      TString tmpStr(objArray->At(0)->GetName());
+      tmpStr.ReplaceAll("p_", "#it{p}_");
+      tmpStr.ReplaceAll("00", "");
+      paveText4->AddText(tmpStr + "GeV/#it{c}");
+
+      TString tmpStr(objArray->At(1)->GetName());
+      tmpStr.ReplaceAll("p_", "#it{p}_");
+      tmpStr.ReplaceAll("00", "");
+      paveText4->AddText(tmpStr + "GeV/#it{c}");
+    
+      paveText4->Draw();
+    
+      //     DrawLatex(0.27, 0.19, 1, Form("%sGeV/#it{c}       %sGeV/#it{c}", objArray->At(0)->GetName(), objArray->At(1)->GetName()), fontSize);
+    
+      gPad->GetCanvas()->SaveAs(Form("ridge_fit_%d_%d.png", i, j));
+      gPad->GetCanvas()->SaveAs(Form("ridge_fit_%d_%d.eps", i, j));
+      gPad->GetCanvas()->SaveAs("fig3b.eps");
+    
+      c3->cd();
+      paveText3 = (TPaveText*) paveText4->Clone();
+      paveText3->Draw();
+      gPad->GetCanvas()->SaveAs(Form("ridge_eta_%d_%d.png", i, j));
+      gPad->GetCanvas()->SaveAs(Form("ridge_eta_%d_%d.eps", i, j));
+      gPad->GetCanvas()->SaveAs("fig3c.eps");
+    }
+  
+  if (gStudySystematic != 50 && gStudySystematic != 51)
+    {
+      Int_t phi1 = 1;
+      Int_t phi2 = proj->FindBin(TMath::Pi()/2-0.001);
+      Int_t phi3 = phi2 + 1;
+      Int_t phi4 = proj->GetNbinsX();
+    }
+  else
+    {
+      Printf(">>>> Using |dphi| < 1.2 <<<<");
+      Int_t phi1 = proj->FindBin(-1.2);;
+      Int_t phi2 = proj->FindBin(1.2);
+      Int_t phi3 = proj->FindBin(TMath::Pi() - 1.2);
+      Int_t phi4 = proj->FindBin(TMath::Pi() + 1.2);
+      Printf("%d %d %d %d", phi1, phi2, phi3, phi4);
+    }
+  
+  Double_t nsYield, asYield, nsYieldE, asYieldE;
+  
+  nsYield = proj->IntegralAndError(phi1, phi2, nsYieldE, "width");
+  nsYield -= min * proj->GetBinWidth(1) * (phi2 - phi1 + 1);
+  
+  asYield = proj->IntegralAndError(phi3, phi4, asYieldE, "width");
+  asYield -= min * proj->GetBinWidth(1) * (phi4 - phi3 + 1);
+  
+  AddPoint(graph[0], xValue1, nsYield, 0, nsYieldE);
+  AddPoint(graph[1], xValue2, asYield, 0, asYieldE);
+  
+  if (nsYieldE > 0 && asYieldE > 0)
+    AddPoint(graph[6], xValue1, nsYield / asYield, 0, nsYield / asYield * TMath::Sqrt(nsYieldE * nsYieldE / nsYield / nsYield + asYieldE * asYieldE / asYield / asYield));
+  
+  AddPoint(graph[9], xValue1, 100.0 * ySubPeak / (etaMax * 2) / nsYield, 0, 0);
+  
+  Printf("Yields: %f +- %f ; %f +- %f", nsYield, nsYieldE, asYield, asYieldE);
+  Printf("Yield ratio: %f   baseline for yield %f",nsYield / asYield,min);
+
+}
+
+Int_t colors[] = { 1, 2, 3, 4, kGray+1, 6, 7, 8 };
+
+void DrawSeveral(Int_t n, const char** graphFiles, Int_t id)
+{
+  ReadGraphs(graphFiles[0]);
+  TGraphErrors*** base = graphs;
+
+  Float_t yMax[] = { 0.1, 0.1, 2, 2, 0.25, 0.25, 4, 1.5, 50, 50, 4, 4, 0.5 };
+  Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
+
+  TString baseName(graphFiles[0]);
+  baseName.ReplaceAll(".", "_");
+  TCanvas* canvas = new TCanvas(Form("%s_%d", baseName.Data(), id), Form("%s_%d", baseName.Data(), id), 800, 600);
+  //   Printf("%p", canvas);
+  gPad->SetGridx();
+  gPad->SetGridy();
+  dummy = new TH2F(Form("hist_%s_%d", graphFiles[0], id), Form(";%s;%s", graphs[0][id]->GetXaxis()->GetTitle(), graphs[0][id]->GetYaxis()->GetTitle()), 100, 0, 60, 100, (yMax[id] > 0) ? 0 : yMax[id], (yMax[id] > 0) ? yMax[id] : 0);
+  dummy->SetStats(0);
+  dummy->DrawCopy();
+  
+  TCanvas* canvas2 = new TCanvas(Form("%s_%d_ratio", baseName.Data(), id), Form("%s_%d", baseName.Data(), id), 800, 600);
+  gPad->SetGridx();
+  gPad->SetGridy();
+  dummy = new TH2F(Form("hist_%s_%d_ratio", graphFiles[0], id), Form(";%s;%s ratio", graphs[0][id]->GetXaxis()->GetTitle(), graphs[0][id]->GetYaxis()->GetTitle()), 100, 0, 60, 100, 0, 2);
+  dummy->SetStats(0);
+  dummy->DrawCopy();
+
+  legend = new TLegend(0.4, 0.6, 0.99, 0.99);
+  legend->SetFillColor(0);
+  
+  for (Int_t fc = 0; fc<n; fc++)
+    {
+      ReadGraphs(graphFiles[fc]);
+    
+      for (Int_t i=0; i<NGraphs; i++)
+       {
+         if (TString(graphFiles[0]).Contains("cms") && i != 2 && i != 5) continue;
+    
+         canvas->cd();
+         //       Printf("%p", canvas);
+         graphs[i][id]->SetMarkerStyle(markers[i%8]);
+         graphs[i][id]->SetMarkerColor(colors[fc]);
+         graphs[i][id]->SetLineColor(colors[fc]);
+         GraphShiftX((TGraphErrors*) graphs[i][id]->DrawClone("PSAME"), 0.5 / n * fc);
+
+         if (fc == 0 && graphs[i][id]->GetN() > 0)
+           legend->AddEntry(graphs[i][id], graphs[i][id]->GetTitle(), "P");
+      
+         if (fc > 0)
+           {
+             canvas2->cd();
+             DivideGraphs(graphs[i][id], base[i][id]);
+             GraphShiftX((TGraphErrors*) graphs[i][id]->DrawClone("PSAME"), 0.5 / n * fc);
+           }
+       }
+    }
+
+  canvas->cd();
+  legend->Draw();
+  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+
+  canvas2->cd();
+  legend->Draw();
+  canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
+}
+
+void CMSRidge()
+{
+  CreateGraphStructure();
+  
+  // scaling done for zyam level difference (note that the values in fig2 have to be multiplied by 2 due to the |dphi| instead of dphi
+  
+  // Fig3a
+  // 0.334516912727  0.0104728088208
+  // 0.721360493366  0.0263651653896
+  // 1.1900579331  0.0301999626238
+  // 1.68965302436  0.024186548724
+  // 2.19286114745  0.0123441101352
+  // 2.69021366723  0.00994871155963
+  // 3.38790257273  0.00434394401877
+
+  //   AddPoint(graphs[0][0], 55.3, 0.0263651653896, 0, 0);
+  AddPoint(graphs[2][0], 55.383968, 0.0301999626238 * 0.72 * 2 / 1.336, 0, 0);
+  AddPoint(graphs[5][0], 55.383968, 0.0123441101352 * 0.154 * 2 / 0.229, 0, 0);
+  
+  // Fig3b
+  // 128.344370861  0.0254622516556
+  // 98.1456953642  0.0154940397351
+  // 56.6225165563  0.00493774834437
+  // 17.6158940397  0.00037880794702
+  // Ref multiplicity for centrality 0.000000 to 3.000000: 55.284322
+  // Ref multiplicity for centrality 3.000000 to 10.000000: 41.124550
+  // Ref multiplicity for centrality 10.000000 to 50.000000: 24.334303
+  // Ref multiplicity for centrality 50.000000 to 100.000000: 8.008294
+
+  //   AddPoint(graphs[2][0], 55.3, 0.0254622516556, 0, 0);
+  AddPoint(graphs[2][0], 40.649288, 0.0154940397351 * 0.5 * 2 / 0.985, 0, 0);
+  AddPoint(graphs[2][0], 23.783224, 0.00493774834437 * 0.29 * 2 / 0.506, 0, 0);
+
+  WriteGraphs("graphs_cms.root");
+}
+
+void ALICECMSMethod()
+{
+  CreateGraphStructure();
+  
+  AddPoint(graphs[2][0], 10, 0.027, 0, 0);
+  AddPoint(graphs[4][0], 10, 0.115 / 1.34 * 1, 0, 0);
+  AddPoint(graphs[5][0], 10, 0.019 / 0.152 * 0.126, 0, 0);
+  
+  AddPoint(graphs[2][4], 10, 0.09, 0, 0);
+  //   AddPoint(graphs[4][4], 10, 0.115, 0, 0);
+  AddPoint(graphs[5][4], 10, 0.137, 0, 0);
+
+  AddPoint(graphs[2][5], 10, 0.04, 0, 0);
+  //   AddPoint(graphs[4][5], 10, 0.046, 0, 0);
+  AddPoint(graphs[5][5], 10, 0.06, 0, 0);
+
+  WriteGraphs("graphs_alice_cmsmethod.root");
+}
+
+void GetSystematic(const char* baseFile = "dphi_corr_pA_121119_3.root",TString fileTag = "graphs_121119",TString fileTagHadrons = "graphs_121119")
+{
+  TString baseFile2="0";
+  fileProj = TFile::Open("dphi_proj.root", "RECREATE");
+  fileProj->Close();
+  
+  gROOT->SetBatch(kTRUE);
+  
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + ".root",fileTagHadrons + ".root");
+  
+  gStudySystematic = 10;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion05.root",fileTagHadrons + "_exclusion05.root");
+  
+  gStudySystematic = 11;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion00.root",fileTagHadrons + "_exclusion00.root");
+
+  gStudySystematic = 12;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusionAS.root",fileTagHadrons + "_exclusionAS.root");
+
+  gStudySystematic = 13;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusionScale.root",fileTagHadrons + "_exclusionScale.root");
+
+  gStudySystematic = 14;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_exclusion12.root",fileTagHadrons + "_exclusion12.root");
+    
+  gStudySystematic = 20;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_nonclosure.root",fileTagHadrons + "_nonclosure.root");
+  
+  gStudySystematic = 30;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_baseline.root",fileTagHadrons + "_baseline.root");
+  
+  gStudySystematic = 40;
+  baseFile2=Form("~/Dropbox/pid-ridge/pPb/dphi_corr_LHC13bc_AOD126_20130417_%s_SystTrkCuts.root",gParticleSpecieName[gParticleSpecie]);
+  CorrelationSubtractionHistogram(baseFile2.Data(), kTRUE, fileTag + "_trackcuts.root",fileTagHadrons + "_trackcuts.root");
+  
+  gStudySystematic = 41;
+  baseFile2=Form("~/Dropbox/pid-ridge/pPb/dphi_corr_LHC13bc_AOD126_20130417_%s_PairCuts.root",gParticleSpecieName[gParticleSpecie]);
+  CorrelationSubtractionHistogram(baseFile2.Data(), kTRUE, fileTag + "_paircuts.root",fileTagHadrons + "_paircuts.root");
+
+  gStudySystematic = 60;
+  CorrelationSubtractionHistogram(baseFile, kTRUE, fileTag + "_otherperipheral.root",fileTagHadrons + "_otherperipheral.root");
+  
+  if(gParticleSpecie!=0)
+    {//PID
+      gStudySystematic = 70;
+      baseFile2=Form("~/Dropbox/pid-ridge/pPb/dphi_corr_LHC13bc_AOD126_20130417_%s_Syst2Sig.root",gParticleSpecieName[gParticleSpecie]);
+      CorrelationSubtractionHistogram(baseFile2.Data(), kTRUE, fileTag + "_syst2sig.root",fileTagHadrons + "_syst2sig.root");
+      
+      gStudySystematic = 71;
+      baseFile2=Form("~/Dropbox/pid-ridge/pPb/dphi_corr_LHC13bc_AOD126_20130417_%s_Syst1Sig.root",gParticleSpecieName[gParticleSpecie]);
+      CorrelationSubtractionHistogram(baseFile2.Data(), kTRUE, fileTag + "_syst1sig.root",fileTagHadrons + "_syst1sig.root");
+      
+      gStudySystematic = 72;
+      baseFile2=Form("~/Dropbox/pid-ridge/pPb/dphi_corr_LHC13bc_AOD126_20130417_%s_SystExclusive.root",gParticleSpecieName[gParticleSpecie]);
+      CorrelationSubtractionHistogram(baseFile2.Data(), kTRUE, fileTag + "_systexclusive.root",fileTagHadrons + "_systexclusive.root");
+      
+    }
+}
+  
+void DrawSystematics(TString fileTag = "graphs_121119")
+{
+  if (0)
+    {
+      const Int_t n = 6;
+      const char* filesBase[] = { "", "_exclusion05", "_exclusion12", "_exclusion00", "_exclusionAS", "_exclusionScale" };
+    }
+  else if(1)
+    {
+      const Int_t n = 6;
+      const char* filesBase[] = { "", "_paircuts" ,"_trackcuts", "_nonclosure", "_baseline", "_otherperipheral" };
+    }
+  else
+    {
+      const Int_t n = 4;
+      const char* filesBase[] = { "", "_syst2sig" ,"_syst1sig", "_systexclusive"};
+    }
+  
+  const char* files[n];
+  for (Int_t i=0; i<n; i++)
+    {
+      str = new TString;
+      str->Form("%s%s.root", fileTag.Data(), filesBase[i]);
+      files[i] = str->Data();
+    }
+  
+  const Int_t plots = 3;
+  Int_t ids[7] = { 0, 1, 4 };
+  
+  for (Int_t i=0; i<plots; i++)
+    {
+      DrawSeveral(n, files, ids[i]);
+      //     break;
+    }
+  
+  new TCanvas;
+  for (Int_t i=0; i<n; i++)
+    DrawLatex(0.2, 0.9 - 0.1 * i, colors[i], files[i]);
+}
+
+void DrawSeveral(const char* graphFile1, const char* graphFile2, const char* graphFile3, Int_t id)
+{
+  const char* files[3] = { graphFile1, graphFile2, graphFile3 };
+  Int_t n = 1;
+  if (graphFile2)
+    n++;
+  if (graphFile3)
+    n++;
+  
+  DrawSeveral(n, files, id);
+}
+  
+void DrawYield(const char* graphFile)
+{
+  DrawGraph(graphFile, 0, 1, "Ridge yield per #Delta#eta", "fig4b.eps", kTRUE);
+}
+
+void DrawRMS(const char* graphFile)
+{
+  DrawGraph(graphFile, 2, 3, "#sigma", 0);
+}
+
+void Drawv2v3(const char* graphFile)
+{
+  DrawGraph(graphFile, 4, 5, "#it{v}_{2} , #it{v}_{3}", "fig4a.eps");
+}
+
+void AddSystUnc(TGraphErrors* graph, Float_t syst020, Float_t syst2060)
+{
+  for (Int_t j=0; j<graph->GetN(); j++)
+    {
+      Float_t syst = syst2060;
+      if (graph->GetX()[j] < 20)
+       syst = syst020;
+  
+      //     Printf("%d %f", j, syst);
+      graph->GetEY()[j] = TMath::Sqrt(graph->GetEY()[j] * graph->GetEY()[j] + syst * syst * graph->GetY()[j] * graph->GetY()[j]);
+    }
+}
+
+void SetSystUnc(TGraphErrors* graph, Float_t syst020, Float_t syst2060)
+{
+  for (Int_t j=0; j<graph->GetN(); j++)
+    {
+      Float_t syst = syst2060;
+      if (graph->GetX()[j] < 20)
+       syst = syst020;
+  
+      //     Printf("%d %f", j, syst);
+      graph->GetEY()[j] = syst * graph->GetY()[j];
+    }
+}
+
+void SetXError(TGraphErrors* graph, Float_t value)
+{
+  for (Int_t j=0; j<graph->GetN(); j++)
+    graph->GetEX()[j] = value;
+}
+
+void DrawGraph(const char* graphFile, Int_t id1, Int_t id2, const char* yLabel = 0, const char* outputFileName, Bool_t corrGraph = kFALSE)
+{
+  ReadGraphs(graphFile);
+  TGraphErrors*** graphsSyst = graphs;
+  ReadGraphs(graphFile);
+  TGraphErrors*** graphsCombinedError = graphs;
+  ReadGraphs(graphFile);
+  
+  if (yLabel == 0)
+    yLabel = graphs[0][id1]->GetYaxis()->GetTitle();
+  
+  Float_t yMax[] = { 0.12, 0.12, 2, 2, 0.2, 0.25, 4, 1.5, 50, 50, 4 };
+  Int_t markers[] = { 20, 21, 22, 23, 29, 33 };
+  Int_t markers2[] = { 24, 25, 26, 32, 30, 27 };
+  Float_t markerSize[] = { 1.7, 1.7, 1.7, 1.7, 2.0, 2.0 };
+
+  if (1)
+    {
+      // default systs (for all pT, centr)
+      Float_t syst020Array[] = { 0.16, 0.18, 0.15, 0.15, 0.14, 0.23 };
+      Float_t syst2060Array[] = { 0.23, 0.25, 0.23, 0.23, 0.18, 0.42 };
+    
+      // 0.5<1 0.5<1
+      Float_t syst020ArrayGr0[] = { 0.23, 0.25, 0.15, 0.15, 0.14, 0.23 };
+      Float_t syst2060ArrayGr0[] = { 0.42, 0.43, 0.23, 0.23, 0.18, 0.42 };
+
+      for (Int_t i=0; i<NGraphs; i++)
+       {
+         Float_t syst020 = syst020Array[id1];
+         Float_t syst2060 = syst2060Array[id1];
+      
+         if (i == 0)
+           {
+             syst020 = syst020ArrayGr0[id1];
+             syst2060 = syst2060ArrayGr0[id1];
+           }
+      
+         //       Printf("%d %d %f %f", i, id1, syst020, syst2060);
+         SetSystUnc(graphsSyst[i][id1], syst020, syst2060);
+         AddSystUnc(graphsCombinedError[i][id1], syst020, syst2060);
+         SetXError(graphsSyst[i][id1], 0.3);
+
+         Float_t syst020 = syst020Array[id2];
+         Float_t syst2060 = syst2060Array[id2];
+      
+         if (i == 0)
+           {
+             syst020 = syst020ArrayGr0[id2];
+             syst2060 = syst2060ArrayGr0[id2];
+           }
+
+         SetSystUnc(graphsSyst[i][id2], syst020, syst2060);
+         AddSystUnc(graphsCombinedError[i][id2], syst020, syst2060);
+         SetXError(graphsSyst[i][id2], 0.3);
+       }
+    
+      //     graphs[0][id1]->Print();
+    }
+  
+  const char* eventClass[] = { "0-20%", "20-40%", "40-60%" };
+  Printf("\n%s:", graphTitles[id1]);
+  for (Int_t i=0; i<3; i++)
+    {
+      Printf("Event class: %s minus 60-100%%", eventClass[i]);
+    
+      for (Int_t j=0; j<NGraphs; j++)
+       if (graphs[j][id1]->GetN() > i)
+         Printf("%s: %.4f +- %.4f (stat) +- %.4f (syst)", graphs[j][id1]->GetTitle(), graphs[j][id1]->GetY()[i], graphs[j][id1]->GetEY()[i], graphsSyst[j][id1]->GetEY()[i]);
+    }
+    
+  Printf("\n%s:", graphTitles[id2]);
+  for (Int_t i=0; i<3; i++)
+    {
+      Printf("Event class: %s minus 60-100%%", eventClass[i]);
+    
+      for (Int_t j=0; j<NGraphs; j++)
+       if (graphs[j][id2]->GetN() > i)
+         Printf("%s: %.4f +- %.4f (stat) +- %.4f (syst)", graphs[j][id2]->GetTitle(), graphs[j][id2]->GetY()[i], graphs[j][id2]->GetEY()[i], graphsSyst[j][id2]->GetEY()[i]);
+    }
+
+  TCanvas* canvas = new TCanvas;
+  gPad->SetTopMargin(0.03);
+  gPad->SetRightMargin(0.01);
+  gPad->SetLeftMargin(0.14);
+  gPad->SetBottomMargin(0.13);
+  //   gPad->SetGridx();
+  //   gPad->SetGridy();
+  TH2F* dummy = new TH2F("dummy", Form(";%s;%s", graphs[0][id1]->GetXaxis()->GetTitle(), yLabel), 3, 0, 60, 100, 0, yMax[id1]);
+  dummy->GetYaxis()->SetTitleOffset(1.1);
+  dummy->SetStats(0);
+  dummy->GetYaxis()->SetNdivisions(505);
+  dummy->GetXaxis()->SetLabelSize(0.06);
+  dummy->GetYaxis()->SetLabelSize(0.06);
+  dummy->GetXaxis()->SetTitleSize(0.06);
+  dummy->GetYaxis()->SetTitleSize(0.06);
+  dummy->Draw();
+  
+  if (strcmp(graphs[0][id1]->GetXaxis()->GetTitle(), "Event class") == 0)
+    {
+      dummy->GetXaxis()->SetBinLabel(1, "0-20%");
+      dummy->GetXaxis()->SetBinLabel(2, "20-40%");
+      dummy->GetXaxis()->SetBinLabel(3, "40-60%");
+      dummy->GetXaxis()->SetLabelSize(0.09);
+      dummy->GetXaxis()->SetTitleOffset(1);
+    }
+  
+  legend = new TLegend((id1 == 4) ? 0.47 : 0.33, (id1 == 4) ? 0.70 : 0.55, 0.95, 0.92);
+  legend->SetNColumns(2);
+  if (id1 == 0 || id1 == 2)
+    legend->SetHeader("Near side   Away side");
+  else if (id1 == 4)
+    legend->SetHeader("  #it{v}_{2}       #it{v}_{3}");
+    
+  legend->SetFillColor(0);
+  legend->SetBorderSize(0);
+
+  if (1)
+    {
+      Int_t fillStyle[11] = { 3001, 3002, 3003, 3004, 3005, 3006, 3007, 3008, 3009, 3010, 3011 };
+    
+      for (Int_t i=0; i<NGraphs; i++)
+       {
+         graphsSyst[i][id1]->SetMarkerStyle(1);
+         graphsSyst[i][id1]->SetMarkerColor(1);
+         graphsSyst[i][id1]->SetLineColor(1);
+         graphsSyst[i][id1]->SetFillColor(1);
+         graphsSyst[i][id1]->SetFillStyle(3001);
+         graphsSyst[i][id1]->Draw("2SAME");
+
+         graphsSyst[i][id2]->SetMarkerStyle(1);
+         graphsSyst[i][id2]->SetMarkerColor(2);
+         graphsSyst[i][id2]->SetLineColor(2);
+         graphsSyst[i][id2]->SetFillColor(2);
+         graphsSyst[i][id2]->SetFillStyle(3001);
+         graphsSyst[i][id2]->Draw("2SAME");
+       }
+    }
+  else
+    Printf(">>>>>>>>>>>> SKIPPING SYST");
+
+  
+  for (Int_t i=0; i<NGraphs; i++)
+    {
+      graphs[i][id1]->SetMarkerStyle(markers[i]);
+      graphs[i][id1]->SetMarkerSize(markerSize[i]);
+      graphs[i][id1]->SetMarkerColor(1);
+      graphs[i][id1]->SetLineColor(1);
+      graphs[i][id1]->Draw("PSAME");
+
+      graphs[i][id2]->SetMarkerStyle(markers2[i]);
+      graphs[i][id2]->SetMarkerSize(markerSize[i]);
+      graphs[i][id2]->SetMarkerColor(2);
+      graphs[i][id2]->SetLineColor(2);
+      graphs[i][id2]->Draw("PSAME");
+
+      TString label(graphs[i][id1]->GetTitle());
+      label.ReplaceAll(".00", ".0");
+      label.ReplaceAll(".50", ".5");
+    
+    
+      /*    label.ReplaceAll(".0", " GeV/#it{c}");
+           label.ReplaceAll(".5", ".5 GeV/#it{c}");*/
+      TObjArray* objArray = label.Tokenize("-");
+      label.Form("%s-%s", objArray->At(0)->GetName(), objArray->At(1)->GetName());
+      label.ReplaceAll("-", ";");
+    
+      if (id1 == 4)
+       {
+         // reduce label
+         label.ReplaceAll(" ", "");
+         label.ReplaceAll(";", "<");
+         tokens = label.Tokenize("<");
+         label.Form("%s < %s < %s < %s ", tokens->At(0)->GetName(), tokens->At(1)->GetName(), tokens->At(4)->GetName(), tokens->At(5)->GetName());
+       }
+    
+      label.ReplaceAll("p_", "#it{p}_");
+      label += "GeV/#it{c}";
+
+      if (graphs[i][id1]->GetN() > 0)
+       {
+         legend->AddEntry(graphs[i][id1], (id1 == 4) ? " " : "    ", "P");
+         legend->AddEntry(graphs[i][id2], label, "P");
+       }
+    }
+  
+  legend->Draw();
+  DrawLatex(0.7, 0.92, 1, "p-Pb #sqrt{s_{NN}} = 5.02 TeV", 0.04);
+  
+  canvas->SaveAs(outputFileName);
+  
+  if (!corrGraph)
+    return;
+  
+  corr = new TGraphErrors;
+  for (Int_t i=0; i<NGraphs; i++)
+    {
+      for (Int_t j=0; j<graphs[i][id1]->GetN(); j++)
+       AddPoint(corr, graphsCombinedError[i][id1]->GetY()[j], graphsCombinedError[i][id2]->GetY()[j], graphsCombinedError[i][id1]->GetEY()[j], graphsCombinedError[i][id2]->GetEY()[j]);
+    }
+
+  new TCanvas;
+  gPad->SetTopMargin(0.03);
+  gPad->SetRightMargin(0.02);
+  gPad->SetLeftMargin(0.14);
+  gPad->SetBottomMargin(0.13);
+  TString titleString;
+  titleString.Form(";%s;%s", graphs[0][id1]->GetYaxis()->GetTitle(), graphs[0][id2]->GetYaxis()->GetTitle());
+  titleString.ReplaceAll("NS", "Near-side");
+  titleString.ReplaceAll("AS", "Away-side");
+  titleString.ReplaceAll("Ridge Yield", "ridge yield per #Delta#eta");
+  dummy = new TH2F("dummy2", titleString, 100, 0, 0.095, 100, 1e-5, 0.095);
+  dummy->GetYaxis()->SetTitleOffset(1.1);
+  dummy->SetStats(0);
+  dummy->GetXaxis()->SetNdivisions(505);
+  dummy->GetYaxis()->SetNdivisions(505);
+  dummy->GetXaxis()->SetLabelSize(0.06);
+  dummy->GetYaxis()->SetLabelSize(0.06);
+  dummy->GetXaxis()->SetTitleSize(0.06);
+  dummy->GetYaxis()->SetTitleSize(0.06);
+  dummy->Draw();
+
+  corr->Draw("PSAME");
+  
+  line = new TLine(0, 0, 0.095, 0.095);
+  line->SetLineWidth(2);
+  line->SetLineStyle(2);
+  line->Draw();
+
+  DrawLatex(0.18, 0.92, 1, "p-Pb #sqrt{s_{NN}} = 5.02 TeV", 0.04);
+}
+
+void PaperCorrFuncAll(const char* fileName)
+{
+  Int_t n = 6;
+  Int_t is[] = { 0, 1, 1, 2, 2, 2, 3 };
+  Int_t js[] = { 1, 1, 2, 1, 2, 3, 3 };
+  Int_t centr[] = { 0, 1, 3, 4 };
+
+  for (Int_t i=0; i<n; i++)
+    for (Int_t j=0; j<4; j++)
+      PaperCorrFunc(fileName, is[i], js[i], centr[j]);
+}
+
+void PaperCorrFunc(const char* fileName, Int_t i = 2, Int_t j = 2, Int_t centr = 0)
+{
+  TFile::Open(fileName);
+  
+  TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centr));
+
+  if (!hist1)
+    return 0;
+  // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+  hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+
+  hist1 = (TH2*) hist1->Clone();
+  hist1->Rebin2D(2, 2); hist1->Scale(0.25);
+  //   hist1->Rebin2D(2, 1); hist1->Scale(0.5);
+  
+  hist1->GetYaxis()->SetRangeUser(-1.99, 1.99);
+  hist1->GetXaxis()->SetTitleOffset(1.5);
+  hist1->GetYaxis()->SetTitleOffset(2);
+  hist1->GetZaxis()->SetTitle(kCorrFuncTitle);
+  hist1->SetStats(kFALSE);
+  hist1->GetZaxis()->SetTitleOffset(2);
+  hist1->GetZaxis()->SetNdivisions(504);
+  hist1->SetStats(kFALSE);
+  hist1->GetZaxis()->CenterTitle(kTRUE);
+  hist1->GetYaxis()->CenterTitle(kTRUE);
+  hist1->GetXaxis()->CenterTitle(kTRUE);
+  hist1->GetYaxis()->SetNdivisions(505);
+  hist1->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+
+  c = new TCanvas("c", "c", 600, 600);
+  gPad->SetPad(0, 0, 1, 1);
+  gPad->SetLeftMargin(0.2);
+
+  TString label(hist1->GetTitle());
+  hist1->SetTitle("");
+  label.ReplaceAll(".00", "");
+  label.ReplaceAll(".0", "");
+  //   label.ReplaceAll(".00", " GeV/#it{c}");
+  //   label.ReplaceAll(".0", " GeV/#it{c}");
+  TObjArray* objArray = label.Tokenize("-");
+  TPaveText* paveText = new TPaveText(0.03, 0.86, 0.44, 0.97, "BRNDC");
+  paveText->SetTextSize(0.035);
+  paveText->SetFillColor(0);
+  paveText->SetShadowColor(0);
+  paveText->SetBorderSize(0);
+  paveText->SetFillStyle(0);
+  
+  TString tmpStr(objArray->At(0)->GetName());
+  tmpStr.ReplaceAll("p_", "#it{p}_");
+  paveText->AddText(tmpStr + "GeV/#it{c}");
+
+  TString tmpStr(objArray->At(1)->GetName());
+  tmpStr.ReplaceAll("p_", "#it{p}_");
+  paveText->AddText(tmpStr + "GeV/#it{c}");
+  //   paveText->AddText(objArray->At(1)->GetName());
+
+  TPaveText* paveText2 = new TPaveText(0.63, 0.86, 0.97, 0.97, "BRNDC");
+  paveText2->SetTextSize(0.035);
+  paveText2->SetBorderSize(0);
+  paveText2->SetFillColor(0);
+  paveText2->SetFillStyle(0);
+  paveText2->SetShadowColor(0);
+  if (objArray->GetEntries() == 4)
+    {
+      paveText2->AddText(Form("p-Pb #sqrt{s_{NN}} = 5.02 TeV"));
+      paveText2->AddText(Form("%s-%s", objArray->At(2)->GetName(), objArray->At(3)->GetName()));
+    }
+  else
+    paveText2->AddText(Form("%s #sqrt{s} = 2.76 TeV", objArray->At(2)->GetName()));
+
+  hist1->Draw("SURF1");
+  paveText->Draw();
+  paveText2->Draw();
+  
+  if (i == 2 && j == 2 && centr == 0)
+    c->SaveAs("fig1b.eps");
+  else if (i == 2 && j == 2 && centr == 4)
+    c->SaveAs("fig1a.eps");
+  
+  c->SaveAs(Form("corr_%d_%d_%d.eps", i, j, centr));
+}
+
+void ExtractSystematics(const char* fileName = "dphi_proj.root")
+{
+  Int_t i = 2; Int_t j = 2;
+  //   Int_t i = 0; Int_t j = 1;
+  Int_t centr = 0;
+  
+  Int_t nSyst = 8;
+  Int_t syst[] = { 0, 10, 11, 12, 13, 20, 30, 40 };
+  
+  TFile::Open(fileName);
+  base = (TH1*) gFile->Get(Form("proj_%d_%d_%d_%d", i, j, centr, syst[0]));
+  
+  c = new TCanvas;
+  //   base->Rebin(2); base->Scale(0.5);
+  base->Draw();
+  
+  Float_t baseValue = base->Integral() / base->GetNbinsX();
+  //   (base->GetBinContent(1) + base->GetBinContent(base->GetNbinsX()/2+1)) / 2;
+  
+  c2 = new TCanvas;
+  c3 = new TCanvas;
+  
+  Int_t color = 2;
+
+  for (Int_t n = 1; n<nSyst; n++)
+    {
+      hist = (TH1*) gFile->Get(Form("proj_%d_%d_%d_%d", i, j, centr, syst[n]));
+      if (!hist)
+       continue;
+    
+      //     hist->Rebin(2); hist->Scale(0.5);
+      c->cd();
+      hist->SetLineColor(color++);
+      Float_t baseValue2 = hist->Integral() / hist->GetNbinsX();
+      Printf("%f %f", baseValue, baseValue2);
+      //     hist->Add(new TF1("func", "1", -5, 5), baseValue - baseValue2);
+      hist->Scale(baseValue / baseValue2);
+      hist->DrawCopy("SAME");
+    
+      c2->cd();
+      hist->DrawCopy((n == 1) ? "HIST" : "HIST SAME")->Divide(base);
+
+      c3->cd();
+      hist->Add(base, -1);
+      hist->DrawCopy((n == 1) ? "HIST" : "HIST SAME")->GetYaxis()->SetRangeUser(-0.1, 0.1);
+    }
+}
+
+void GetProjectionSystematics(Float_t eta)
+{
+  fileProj = TFile::Open("dphi_proj2.root", "RECREATE");
+  fileProj->Close();
+
+  Int_t i = 2; Int_t j = 2;
+  //   Int_t i = 0; Int_t j = 1;
+
+  Int_t n = 5;
+  const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121116_global.root", "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121116_hybrid.root" /*pair cuts*/, "dphi_corr_pA_121119_hijing.root" };
+  
+  for (Int_t centr=0; centr<6; centr++)
+    {
+      for (Int_t k=0; k<n; k++)
+       {
+         gStudySystematic = 0;
+         if (k == 2)
+           gStudySystematic = 20;
+      
+         TFile::Open(files[k]);
+         const char* label = 0;
+         TH1* hist = GetProjections(i, j, centr, &label, eta);
+         if (!hist)
+           continue;
+      
+         fileProj = TFile::Open("dphi_proj2.root", "UPDATE");
+         hist->Write(Form("proj_%d_%d_%d_%d", i, j, centr, k*10));
+         fileProj->Close();
+       }
+    }
+}
+
+void CMSPlot()
+{
+  Int_t centrArr[] = { 0, 1, 3, 4 };
+  Int_t nCentr = 4;
+  
+  //   Int_t i = 1; Int_t j = 2;
+  Int_t i = 2; Int_t j = 2;
+
+  if (1)
+    {
+      const char* labels[] = { "Data", "HIJING", "DPMJET" };
+      //   const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_121121_hijing_uncorrected.root", "dphi_corr_121121_dpmjet_uncorrected.root" };
+      //     const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_121121_hijing_step0.root", "dphi_corr_121121_dpmjet_step0.root" };
+      const char* files[] = { "dphi_corr_pA_121121_cmsmethod.root", "dphi_corr_121122_cmsmethod_hijing_step0.root", "dphi_corr_121122_cmsmethod_dpmjet_step0.root" };
+      //     const char* files[] = { "dphi_corr_pA_121122_cnd_cmsmethod.root", "dphi_corr_121122_cmsmethod_hijing_cnd_step0.root", 0 };
+      Int_t nFiles = 3;
+    }
+  else
+    {
+      const char* labels[] = { "ALICE method", "CMS method" };
+      //     const char* files[] = { "dphi_corr_pA_121119_3.root", "dphi_corr_pA_121121_cmsmethod.root" };
+      const char* files[] = { "dphi_corr_121121_hijing_step0.root", "dphi_corr_121122_cmsmethod_hijing_step0.root" };
+      Int_t nFiles = 2;
+    }
+  
+  Int_t colors[] = { 1, 2, 4 };
+  
+  c = new TCanvas("c", "c", 800, 800);
+  c->Divide(2, 2);
+  
+  Float_t maxValue = -1;
+  for (Int_t centr = 0; centr < nCentr; centr++)
+    {
+      c->cd(centr+1);
+    
+      legend = new TLegend(0.15, 0.55, 0.46, 0.85);
+      legend->SetFillColor(0);
+    
+      for (Int_t fileId = 0; fileId < nFiles; fileId++)
+       {
+         if (files[fileId] == 0)
+           continue;
+      
+         TFile::Open(files[fileId]);
+  
+         TH2* hist1 = (TH2*) gFile->Get(Form("dphi_%d_%d_%d", i, j, centrArr[centr]));
+         if (!hist1)
+           return 0;
+         hist1 = (TH2*) hist1->Clone(Form("%s_%.1f", hist1->GetName(), 0.0));
+
+         // NOTE fix normalization. these 2d correlations come out of AliUEHist normalized by dphi bin width, but not deta
+         hist1->Scale(1.0 / hist1->GetYaxis()->GetBinWidth(1));
+      
+         tokens = TString(hist1->GetTitle()).Tokenize("-");
+         centralityStr = new TString;
+         if (tokens->GetEntries() > 2)
+           *centralityStr = tokens->At(2)->GetName();
+         if (tokens->GetEntries() > 3)
+           *centralityStr = *centralityStr + "-" + tokens->At(3)->GetName();
+      
+         Float_t etaMin = 1.0;
+      
+         proj1x = hist1->ProjectionX(Form("proj1x_%d_%d_%d_%d", i, j, centr, fileId), hist1->GetYaxis()->FindBin(-etaMax+0.01), hist1->GetYaxis()->FindBin(-etaMin-0.01));
+         proj2x = hist1->ProjectionX(Form("proj2x_%d_%d_%d_%d", i, j, centr, fileId), hist1->GetYaxis()->FindBin(etaMin+0.01), hist1->GetYaxis()->FindBin(etaMax-0.01));
+         proj1x->Add(proj2x);
+      
+         proj1x->Scale(hist1->GetYaxis()->GetBinWidth(1));
+         proj1x->Scale(1.0 / (etaMax - etaMin) / 2);
+
+         proj1x->Rebin(2); proj1x->Scale(0.5);
+         proj1x->Rebin(2); proj1x->Scale(0.5);
+      
+         // zyam
+         Float_t zyam = proj1x->Integral(proj1x->FindBin(TMath::Pi()/2 - 0.5), proj1x->FindBin(TMath::Pi()/2 + 0.5)) / (proj1x->FindBin(TMath::Pi()/2 + 0.5) - proj1x->FindBin(TMath::Pi()/2 - 0.5) + 1);
+         proj1x->Add(new TF1("f", "1", -5, 5), -zyam);
+
+         proj1x->SetStats(kFALSE);
+         proj1x->GetXaxis()->SetTitleOffset(1);
+      
+         proj1x->SetLineColor(colors[fileId]);
+         if (maxValue < 0)
+           maxValue = proj1x->GetMaximum() * 2;
+         proj1x->SetMaximum(maxValue);
+         //       proj1x->SetMinimum(-0.01);
+         proj1x->Draw((fileId == 0) ? "" : "SAME");
+      
+         legend->AddEntry(proj1x, labels[fileId]);
+       }
+      legend->Draw();
+    }
+}
+
+void FourierFactorization(const char* fileName)
+{
+  ReadGraphs(fileName);
+  
+  Int_t n = 6;
+  Int_t is[] =    { 0, 1, 1, 2, 2, 2, 3 };
+  Int_t js[] =    { 1, 1, 2, 1, 2, 3, 3 };
+  Bool_t symm[] = { 1, 0, 1, 0, 0, 1, 0 };
+
+  Int_t graphID = 5;
+  
+  TGraphErrors** symmGraphs[] = { graphs[0], graphs[2], graphs[5] };
+  
+  for (Int_t i=0; i<symmGraphs[0][graphID]->GetN(); i++)
+    {
+      Printf("%d", i);
+      graphs[1][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[0][graphID]->GetY()[i] * symmGraphs[1][graphID]->GetY()[i]);
+      graphs[3][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[2][graphID]->GetY()[i] * symmGraphs[1][graphID]->GetY()[i]);
+      graphs[4][graphID]->GetY()[i] = TMath::Sqrt(symmGraphs[2][graphID]->GetY()[i] * symmGraphs[2][graphID]->GetY()[i]);
+    }
+  
+  //   graphs[1][graphID]->Draw("A*");
+  
+  WriteGraphs();
+}
+
+void CompareATLAS()
+{
+  atlas = ReadHepdata("/home/jgrosseo/Dropbox/alice-paper-paridge/comparison_atlas/atlas_figaux10.txt", kFALSE, 1, 1);
+  atlas->SetMarkerStyle(20);
+  atlas->Draw("PA");
+  atlas->SetTitle("");
+  atlas->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  atlas->GetYaxis()->SetTitle("s_{2} / v_{2}");
+  atlas->GetYaxis()->SetRangeUser(0.02, 0.18);
+
+  alice = new TGraphErrors;
+  AddPoint(alice, 0.75, 0.0583503, 0.25, 0.00831755);
+  AddPoint(alice, 1.5, 0.0953562, 0.5, 0.0134597);
+  AddPoint(alice, 3, 0.128309, 1, 0.0189634);
+  
+  alice->SetMarkerStyle(21);
+  alice->SetLineColor(2);
+  alice->SetMarkerColor(2);
+  alice->Draw("PSAME");
+}
+
+Double_t CalculateVnAnalyticallyRaw(TH1* proj,Int_t n, Float_t BaseLineShift){
+  //proj->Print("all");
+  Double_t num=0;
+  Double_t den=0;
+  Printf("BaseLineShift %f",BaseLineShift);
+  for(Int_t ibin=1;ibin<=proj->GetNbinsX();ibin++){
+    if(proj->GetBinContent(ibin)==0)continue;
+    num+=TMath::Cos(n*proj->GetBinCenter(ibin))*(proj->GetBinContent(ibin)+BaseLineShift);
+    den+=proj->GetBinContent(ibin)+BaseLineShift;
+    //Printf("num %f   den %f BaselineShift %f",num,den,BaseLineShift);
+  }
+  Double_t vn=num/den;
+  return vn;
+}
+
+Double_t CalculateVnAnalytically(TH1* proj,Int_t n, Float_t BaseLineShift){
+  Double_t vn=CalculateVnAnalyticallyRaw(proj, n, BaseLineShift);
+  if (vn < 0)
+    return 0;
+  Printf("v2 Value: %f",TMath::Sqrt(vn));
+  return TMath::Sqrt(vn);
+}
+
+Double_t CalculateVnErrAnalytically(TH1* proj,Int_t n, Float_t BaseLineShift, Double_t BaseLineShiftError){  
+  Double_t I=0;
+  Double_t v2=CalculateVnAnalytically(proj,n,BaseLineShift);
+  if (v2 <= 0)
+    return 0;
+  
+  for(Int_t ibin=1;ibin<=proj->GetNbinsX();ibin++){
+    if(proj->GetBinContent(ibin)==0)continue;
+    I+=proj->GetBinContent(ibin)+BaseLineShift;
+  }
+  Double_t sigma2=0;
+  for(Int_t ibin=1;ibin<=proj->GetNbinsX();ibin++){
+    if(proj->GetBinContent(ibin)==0)continue;
+    Double_t ck=TMath::Cos(n*proj->GetBinCenter(ibin));
+    sigma2+=(ck-v2)*(ck-v2)*(proj->GetBinError(ibin)+BaseLineShiftError)*(proj->GetBinError(ibin)+BaseLineShiftError);
+  }
+  sigma2=sigma2/(I*I);
+  
+  Double_t vnerr=TMath::Sqrt(sigma2);
+  //we take the sqrt of vn
+  vnerr=vnerr/(2*v2);
+  return vnerr;
+}
+
+
+void UnfoldUsingVnHadrons(TGraphErrors*** Inputgraphs,const char* graphFileHadrons,Int_t id)
+{
+  //create a new set of graphs
+  TGraphErrors*** graphsHadrons = new TGraphErrors**[NGraphs];
+  for (Int_t i=0; i<NGraphs; i++)
+    {
+      graphsHadrons[i] = new TGraphErrors*[NHists];
+      for (Int_t j=0; j<NHists; j++)
+       {
+         graphsHadrons[i][j] = new TGraphErrors;
+         graphsHadrons[i][j]->GetYaxis()->SetTitle(graphTitles[j]);
+       }
+    }
+  //get the hadron ones
+  TFile* fileHadrons = TFile::Open(graphFileHadrons);
+  for (Int_t i=0; i<NGraphs; i++)
+    for (Int_t j=0; j<NHists; j++)
+      graphsHadrons[i][j] = (TGraphErrors*) fileHadrons->Get(Form("graph_%d_%d", i, j));
+  
+  for (Int_t i=0; i<NGraphs; i++){
+    Printf("\nUnfolding %s with hadron file %s",Inputgraphs[i][id]->GetTitle(),graphsHadrons[i][id]->GetTitle());
+    Printf("\033[1;51mInput graph\033[m");
+    Inputgraphs[i][id]->Print("all");
+    Printf("\033[1;51mHadron graph\033[m");
+    graphsHadrons[i][id]->Print("all");
+    for(Int_t ipoint=0;ipoint<Inputgraphs[i][id]->GetN();ipoint++){//loop over the input file
+      Double_t y=0,x=0,xHadrons=0,yHadrons=0;
+      Double_t ey=0,ex=0,exHadrons=0,eyHadrons=0;
+      Inputgraphs[i][id]->GetPoint(ipoint,x,y);
+      Int_t ipointHad=0;
+      Bool_t NoHad=0;
+      graphsHadrons[i][id]->GetPoint(ipointHad,xHadrons,yHadrons);
+      while(xHadrons!=x){//Get the same bin of hadrons, doesn't take directly ipoint because some points might be missing in one of the two histos
+       ipointHad++;
+       graphsHadrons[i][id]->GetPoint(ipointHad,xHadrons,yHadrons);
+       if(ipointHad>10){//the hadron point is missing
+         NoHad=1;
+         break;
+       }
+      }
+      if(NoHad)continue;//if hadron measurement is missing skip the point (unfolding not possible)
+      ex=Inputgraphs[i][id]->GetErrorX(ipoint);
+      ey=Inputgraphs[i][id]->GetErrorY(ipoint);
+      exHadrons=graphsHadrons[i][id]->GetErrorX(ipointHad);
+      eyHadrons=graphsHadrons[i][id]->GetErrorY(ipointHad);
+      printf("(v%d) Input file: %f  %f +- %f  Hadrons: %f  %f +- %f",(id==4)?2:3,x,y,ey,xHadrons,yHadrons,eyHadrons);
+      Double_t unfoldedVal=(y*y)/yHadrons;
+      Inputgraphs[i][id]->SetPoint(ipoint,x,unfoldedVal);
+      Inputgraphs[i][id]->SetPointError(ipoint,ex,unfoldedVal*TMath::Sqrt((2*ey/y)*(2*ey/y)+(eyHadrons/yHadrons)*(eyHadrons/yHadrons)));
+      Inputgraphs[i][id]->GetPoint(ipoint,x,y);
+      ex=Inputgraphs[i][id]->GetErrorX(ipoint);
+      ey=Inputgraphs[i][id]->GetErrorY(ipoint);
+      Printf("  AfterCorrection:  %f  %f +- %f",x,y,ey);
+    }
+  }
+  delete graphsHadrons;
+  fileHadrons->Close();
+  delete fileHadrons;
+}
+
+TString gReferenceFile;
+void DivideOutReferenceParticle(const char* graphFile, const char* outputFile = "graphs.root")
+{
+  ReadGraphs(gReferenceFile);
+  TGraphErrors*** graphsHadrons = graphs;
+  
+  ReadGraphs(graphFile);
+
+  for (Int_t graphID = 4; graphID <= 4; graphID++) //v2 + v3
+  {
+    const Int_t posFullPt = 6; // position where vn for unfolding is located (i.e. full pT,a and pT,t range)
+    Bool_t allPt = TString(graphFile).Contains("allpt");
+    Int_t maxPt = (allPt) ? posFullPt : NGraphs;
+    for (Int_t i=0; i<maxPt; i++)
+    {
+      Int_t referencePos = (allPt) ? posFullPt : i;
+      for (Int_t j=0; j<graphs[i][graphID]->GetN(); j++)
+      {
+       printf("%d %d: %f %f ", i, j, graphs[i][graphID]->GetY()[j], graphsHadrons[referencePos][graphID]->GetY()[j]);
+       if (graphs[i][graphID]->GetY()[j] <= 1e-6)
+       {
+         Printf("");
+         continue;
+       }
+       if (graphsHadrons[referencePos][graphID]->GetY()[j] <= 1e-6)
+       {
+         graphs[i][graphID]->GetY()[j] = 0;
+         graphs[i][graphID]->GetEY()[j] = 0;
+         Printf("");
+         continue;
+       }
+       Float_t relError = graphs[i][graphID]->GetEY()[j] / graphs[i][graphID]->GetY()[j];
+       graphs[i][graphID]->GetY()[j] = graphs[i][graphID]->GetY()[j] * graphs[i][graphID]->GetY()[j] / graphsHadrons[referencePos][graphID]->GetY()[j];
+       graphs[i][graphID]->GetEY()[j] = graphs[i][graphID]->GetY()[j] * TMath::Sqrt(TMath::Power(2.0 * relError, 2) + TMath::Power(graphsHadrons[referencePos][graphID]->GetEY()[j] / graphsHadrons[referencePos][graphID]->GetY()[j], 2));
+       Printf("%f +- %f", graphs[i][graphID]->GetY()[j], graphs[i][graphID]->GetEY()[j]);
+      }
+    }
+  }
+
+  WriteGraphs(outputFile);
+}
+
+void ExtractForSkalarProductComparison(TString inputFileNameBase, TString outputFileNameBase, Bool_t divideOut)
+{
+  gDoSubtraction = 1;
+//   gDoSubtraction = 0; //gUseAnalyticalVn = 1;
+
+  if (1)
+  {
+    gDoEtaGap = 1; 
+//     gDoEtaGapAS = 1;
+    CorrelationSubtractionHistogram(inputFileNameBase + ".root", kTRUE, outputFileNameBase + ".root", 0);
+    if (divideOut)
+      DivideOutReferenceParticle(outputFileNameBase + ".root", outputFileNameBase + "_unfolded.root");
+  }
+  else
+  {
+    gDoEtaGap = 0; 
+  //   gDoEtaGapAS = 0;
+    CorrelationSubtractionHistogram(inputFileNameBase + ".root", kTRUE, outputFileNameBase + "_nogap.root", 0);
+    if (divideOut)
+      DivideOutReferenceParticle(outputFileNameBase + "_nogap.root", outputFileNameBase + "_nogap_unfolded.root");
+  }
+}
+
+void ExtractForSkalarProductComparisonAllSpecies()
+{
+  gROOT->SetBatch(kTRUE);
+
+  TString baseInput("dphi_corr_LHC13bc_20130604");
+  TString baseOutput("graphs_130618_otherbaseline");
+//   baseOutput += "_nosub";
+//   TString postfix("_wingremoved");
+  TString postfix;
+  
+  if (1)
+  {
+    gBinning = 0;
+    gReferenceFile = baseOutput + ".root";
+    ExtractForSkalarProductComparison(baseInput + "_Hadrons" + postfix, baseOutput, kTRUE); 
+    ExtractForSkalarProductComparison(baseInput + "_Pions" + postfix, baseOutput + "_pions", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_Protons" + postfix, baseOutput + "_protons", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_Kaons" + postfix, baseOutput + "_kaons", kTRUE);
+  }
+  
+  if (0)
+  {
+    gBinning = 1;
+    baseInput += "_allpt";
+    baseOutput += "_allpt";
+    gReferenceFile = baseOutput + ".root";
+    ExtractForSkalarProductComparison(baseInput + "_55_Hadrons_symmetrized" + postfix, baseOutput, kTRUE); 
+    ExtractForSkalarProductComparison(baseInput + "_63_Pions" + postfix, baseOutput + "_pions", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_63_Protons" + postfix, baseOutput + "_protons", kTRUE);
+    ExtractForSkalarProductComparison(baseInput + "_63_Kaons" + postfix, baseOutput + "_kaons", kTRUE);
+  }
+}
+
+void Test()
+{
+  TString baseOutput("graphs_130503");
+  const char* arr[] = { "", "_pions", "_protons", "_kaons"};
+  
+  for (Int_t i=1; i<4; i++)
+  {
+    ReadGraphs(baseOutput + arr[i] + ".root");
+    UnfoldUsingVnHadrons(graphs, baseOutput + ".root", 4);
+    WriteGraphs(baseOutput + arr[i] + "_unfolded2.root");
+  }
+}