PWG3/vertexingHF/charmFlow/Extractv2from2Dhistos.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TInterpreter.h>
   3 #include <TString.h>
   4 #include <TObjString.h>
   5 #include <TObjArray.h>
   6 #include <TMath.h>
   7 #include <TFile.h>
   8 #include <TCanvas.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH1D.h>
  12 #include <TF1.h>
  13 #include <TStyle.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TLatex.h>
  17 #include <TPaveText.h>
  18 #include <TDatabasePDG.h>
  19 #include <TGraphErrors.h>
  20 #include <TGraphAsymmErrors.h>
  21
  22
  23 #include "AliAODRecoDecayHF.h"
  24 #include "AliRDHFCuts.h"
  25 #include "AliRDHFCutsDplustoKpipi.h"
  26 #include "AliRDHFCutsD0toKpi.h"
  27 #include "AliHFMassFitter.h"
  28 #include "AliEventPlaneResolution.h"
  29 #endif
  30
  31
  32 // Common variables: to be configured by the user
  33 TString filname="AnalysisResultsgoodruns.root";
  34 TString listname="coutputv2D0Std";//"coutputv2D0RAACuts"; //"coutputv2D0Std";
  35 TString filcutsname="AnalysisResultsgoodruns.root";
  36 Int_t minCent=30;
  37 Int_t maxCent=50;
  38 Int_t mesonPDG=421;
  39
  40 const Int_t nFinalPtBins=3;
  41 Double_t ptlims[nFinalPtBins+1]={2.,5.,8.,12.};
  42 Double_t sigmaRangeForSig=2.5;
  43 Double_t sigmaRangeForBkg=4.5;
  44 Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2};
  45 Bool_t useConstantvsBkgVsMass=kFALSE;
  46 Int_t rebinHistov2Mass[nFinalPtBins]={2,2,5};
  47 Int_t factor4refl=0;
  48 Int_t minPtBin[nFinalPtBins]={-1,-1,-1};
  49 Int_t maxPtBin[nFinalPtBins]={-1,-1,-1};
  50
  51 // systematic errors for 2-5, 5-8 and 8-12 (no run-by-run weights)
  52 /*
  53 Double_t systErrMeth1[nFinalPtBins]={
  54   (0.308-0.169)/2.,
  55   (0.14-0.1)/2.,
  56   (0.04+0.02)/2.
  57 };
  58 Double_t systErrMeth2[nFinalPtBins]={
  59   (0.305-0.252)/2.,
  60   (0.129-0.020)/2.,
  61   (0.101+0.06)/2.
  62 };
  63 */
  64
  65 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights)
  66 Double_t systErrMeth1[nFinalPtBins]={
  67   (0.23-0.10)/2.,
  68   (0.152-0.078)/2.,
  69   (0.161-0.097)/2.
  70 };
  71 Double_t systErrMeth2[nFinalPtBins]={
  72   (0.164-0.097)/2.,
  73   (0.110-0.012)/2.,
  74   (0.131-0.036)/2.
  75 };
  76
  77
  78 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights, RAA cuts)
  79 /*
  80 Double_t systErrMeth1[nFinalPtBins]={
  81   (0.265-0.122)/2.,
  82   (0.165-0.117)/2.,
  83   (0.238-0.169)/2.
  84 };
  85 Double_t systErrMeth2[nFinalPtBins]={
  86   (0.174-0.135)/2.,
  87   (0.18-0.11)/2.,
  88   (0.311-0.28)/2.
  89 };
  90 */
  91
  92 // output of mass fitter
  93 Int_t hMinBin;
  94 Int_t hMaxBin;
  95 Double_t signalFromFit,esignalFromFit;
  96 Double_t massFromFit,sigmaFromFit;
  97 TF1* fB1=0x0;
  98 TF1* fB2=0x0;
  99 TF1* fSB=0x0;
 100
 101 void LoadMassHistos(TList* lst, TH2F** hMassDphi);
 102 Bool_t DefinePtBins(TDirectoryFile* df);
 103 Double_t v2vsMass(Double_t *x, Double_t *par);
 104 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh);
 105 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad);
 106 TH1F* DoSideBands(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c1, Int_t optBkg=0);
 107 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors=0, Bool_t useConst=kTRUE);
 108
 109
 110 void Extractv2from2Dhistos(){
 111
 112   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 113
 114   TFile* filcuts=new TFile(filcutsname.Data());
 115   TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
 116   Bool_t binOK=DefinePtBins(dfcuts);
 117   if(!binOK) return;
 118
 119   TFile* fil=new TFile(filname.Data());
 120   TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
 121   TList* lst=(TList*)df->Get(listname.Data());
 122   if(!lst){
 123     printf("ERROR: list %s not found in file\n",listname.Data());
 124     return;
 125   }
 126   Double_t rcfmin,rcfmax;
 127   Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
 128   Double_t resolSyst=(rcfmax-rcfmin)/2./resolFull;
 129   printf("Relative Systematic Error on RCF=%f\n",resolSyst);
 130
 131   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
 132   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 133     hMassDphi[iFinalPtBin]=0x0;
 134   }
 135   LoadMassHistos(lst, hMassDphi);
 136
 137   TCanvas** c1=new TCanvas*[nFinalPtBins];
 138   TCanvas** c2=new TCanvas*[nFinalPtBins];
 139
 140   Double_t v2M1[nFinalPtBins],errv2M1[nFinalPtBins];
 141   Double_t v2M2[nFinalPtBins],errv2M2[nFinalPtBins];
 142
 143   gStyle->SetPalette(1);
 144   gStyle->SetOptTitle(0);
 145
 146   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 147     printf("**************** BIN %d ******************\n",iFinalPtBin);
 148     printf("\n--------- Method 1: Side Bands ----------\n");
 149     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 150
 151     c1[iFinalPtBin]=new TCanvas(Form("cMeth1PtBin%d",iFinalPtBin),Form("Meth1PtBin%d",iFinalPtBin));
 152     c1[iFinalPtBin]->Divide(2,2);
 153     c1[iFinalPtBin]->cd(1);
 154     hMassDphi[iFinalPtBin]->Draw("colz");
 155     c1[iFinalPtBin]->cd(2);
 156     Bool_t out=FitMassSpectrum(hMass,(TPad*)gPad);
 157     if(!out) continue;
 158
 159     TH1F* hCos2PhiSig=DoSideBands(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],c1[iFinalPtBin]);
 160     Double_t v2obsM1=hCos2PhiSig->GetMean();
 161     Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 162     printf("v2obs = %f +- %f\n",v2obsM1,errv2obsM1);
 163     v2M1[iFinalPtBin]=v2obsM1/resolFull;
 164     errv2M1[iFinalPtBin]=errv2obsM1/resolFull;
 165     printf("v2 = %f +- %f\n",v2M1[iFinalPtBin],errv2M1[iFinalPtBin]);
 166
 167     printf("\n--------- Method 2: S/S+B ----------\n");
 168     c2[iFinalPtBin]=new TCanvas(Form("cMeth2Bin%d",iFinalPtBin),Form("Meth2Bin%d",iFinalPtBin));
 169     TF1* fv2=DoFitv2vsMass(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],c2[iFinalPtBin],0,useConstantvsBkgVsMass);
 170
 171     Double_t v2obsM2=fv2->GetParameter(3);
 172     Double_t errv2obsM2=fv2->GetParError(3);
 173     printf("v2obs = %f +- %f\n",v2obsM2,errv2obsM2);
 174     v2M2[iFinalPtBin]=v2obsM2/resolFull;
 175     errv2M2[iFinalPtBin]=errv2obsM2/resolFull;
 176     printf("v2 = %f +- %f\n",v2M2[iFinalPtBin],errv2M2[iFinalPtBin]);
 177     c1[iFinalPtBin]->SaveAs(Form("cMeth1Bin%d.root",iFinalPtBin));
 178     c2[iFinalPtBin]->SaveAs(Form("cMeth2Bin%d.root",iFinalPtBin));
 179   }
 180
 181   printf("\n--------- Summary ------------\n");
 182
 183   TH1F* hv2m1=new TH1F("hv2m1","Side Band subtraction",nFinalPtBins,ptlims);
 184   TH1F* hv2m2=new TH1F("hv2m2","Fit to v2 vs. mass",nFinalPtBins,ptlims);
 185    for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 186     printf("PtBin %d   v2method1 = %f +- %f   v2method2 = %f +-%f\n",iFinalPtBin,
 187            v2M1[iFinalPtBin],errv2M1[iFinalPtBin],
 188            v2M2[iFinalPtBin],errv2M2[iFinalPtBin]
 189            );
 190     hv2m1->SetBinContent(iFinalPtBin+1,v2M1[iFinalPtBin]);
 191     hv2m1->SetBinError(iFinalPtBin+1,errv2M1[iFinalPtBin]);
 192     hv2m2->SetBinContent(iFinalPtBin+1,v2M2[iFinalPtBin]);
 193     hv2m2->SetBinError(iFinalPtBin+1,errv2M2[iFinalPtBin]);
 194   }
 195
 196    TH1F* hSystErr1=(TH1F*)hv2m1->Clone("hSystErr1");
 197    TH1F* hSystErr2=(TH1F*)hv2m2->Clone("hSystErr2");
 198    for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 199      Double_t systRes1=resolSyst*v2M1[iFinalPtBin];
 200      Double_t systRes2=resolSyst*v2M2[iFinalPtBin];
 201      printf("%f\n",systRes1);
 202      Double_t syste1=TMath::Sqrt(systErrMeth1[iFinalPtBin]*systErrMeth1[iFinalPtBin]+systRes1*systRes1);
 203      Double_t syste2=TMath::Sqrt(systErrMeth2[iFinalPtBin]*systErrMeth2[iFinalPtBin]+systRes2*systRes2);
 204      hSystErr1->SetBinError(iFinalPtBin+1,syste1);
 205      hSystErr2->SetBinError(iFinalPtBin+1,syste2);
 206    }
 207
 208
 209    Double_t maxy=TMath::Max(hv2m2->GetMaximum(),hv2m1->GetMaximum())+0.1;
 210    Double_t miny=TMath::Min(hv2m2->GetMinimum(),hv2m1->GetMinimum())-0.1;
 211    TH2F* hempty=new TH2F("hempty","",10,0.,hv2m1->GetXaxis()->GetXmax()+2.,10,miny,maxy);
 212    hempty->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 213    hempty->GetYaxis()->SetTitle("v_{2}");
 214
 215    TCanvas* cv2=new TCanvas("cv2","v2");
 216    hempty->Draw();
 217    hv2m1->SetMarkerStyle(26);
 218    hSystErr1->SetFillColor(kGray);
 219    hSystErr1->SetFillStyle(3002);
 220    hSystErr1->Draw("e2same");
 221
 222    hv2m2->SetLineColor(kRed+1);
 223    hv2m2->SetMarkerColor(kRed+1);
 224    hv2m2->SetMarkerStyle(20);
 225    hSystErr2->SetFillColor(kRed-9);
 226    hSystErr2->SetFillStyle(3005);
 227    hSystErr2->Draw("e2same");
 228
 229    hv2m1->Draw("same");
 230    hv2m2->Draw("same");
 231
 232    TLegend* leg2=new TLegend(0.5,0.7,0.89,0.89);
 233    leg2->SetFillStyle(0);
 234    TLegendEntry* ent=leg2->AddEntry(hv2m1,"Side Band subtraction","P");
 235    ent->SetTextColor(hv2m1->GetMarkerColor());
 236    ent=leg2->AddEntry(hv2m2,"Fit to v2 vs. mass","P");
 237    ent->SetTextColor(hv2m2->GetMarkerColor());
 238    leg2->Draw();
 239    cv2->Update();
 240    cv2->SaveAs("Dzero-v2-2Dmethods.gif");
 241
 242    TFile* outfil=new TFile("Dzero-v2-2Dmethods.root","recreate");
 243    outfil->cd();
 244    hv2m1->Write();
 245    hv2m2->Write();
 246    hSystErr1->Write();
 247    hSystErr2->Write();
 248    outfil->Close();
 249 }
 250
 251
 252 Double_t v2vsMass(Double_t *x, Double_t *par){
 253   // Fit function for signal+background
 254   // par[0] = S/B at the mass peak
 255   // par[1] = mass
 256   // par[2] = sigma
 257   // par[3] = v2sig
 258   // par[4] = v2back
 259
 260   Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
 261   Double_t fracbkg=1-fracsig;
 262   return par[3]*fracsig+par[4]*fracbkg;
 263 }
 264
 265 Double_t v2vsMassLin(Double_t *x, Double_t *par){
 266   // Fit function for signal+background
 267   // par[0] = S/B at the mass peak
 268   // par[1] = mass
 269   // par[2] = sigma
 270   // par[3] = v2sig
 271   // par[4] = v2back constant
 272   // par[5] = v2back slope
 273
 274   Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
 275   Double_t fracbkg=1-fracsig;
 276   return par[3]*fracsig+(par[4]+par[5]*x[0])*fracbkg;
 277 }
 278
 279 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
 280
 281   Int_t nMassBins=hMass->GetNbinsX();
 282   hMinBin=3;
 283   hMaxBin=nMassBins-2;
 284   Double_t hmin=hMass->GetBinLowEdge(hMinBin);
 285   Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
 286   Float_t massD=TDatabasePDG::Instance()->GetParticle(mesonPDG)->Mass();
 287
 288   AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
 289   fitter->SetReflectionSigmaFactor(factor4refl);
 290   fitter->SetInitialGaussianMean(massD);
 291   Bool_t out=fitter->MassFitter(0);
 292   if(!out) return kFALSE;
 293   fitter->Signal(sigmaRangeForSig, signalFromFit,esignalFromFit);
 294   massFromFit=fitter->GetMean();
 295   sigmaFromFit=fitter->GetSigma();
 296   fB1=fitter->GetBackgroundFullRangeFunc();
 297   fB2=fitter->GetBackgroundRecalcFunc();
 298   fSB=fitter->GetMassFunc();
 299   if(!fB1) return kFALSE;
 300   if(!fB2) return kFALSE;
 301   if(!fSB) return kFALSE;
 302   if(pad){
 303     fitter->DrawHere(gPad,3.,0.);
 304   }
 305   return kTRUE;
 306 }
 307
 308 TH1F* DoSideBands(Int_t iFinalPtBin,
 309                   TH2F* hMassDphi,
 310                   TH1F* hMass,
 311                   Int_t rebin,
 312                   TCanvas* c1,
 313                   Int_t optBkg){
 314
 315   // Build histo with cos(2*deltaphi) distribution for signal
 316
 317   Double_t hmin=hMass->GetBinLowEdge(hMinBin);
 318   Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
 319
 320   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 321   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 322   Int_t minBinSig=hMass->FindBin(minMassSig);
 323   Int_t maxBinSig=hMass->FindBin(maxMassSig);
 324   Double_t minMassSigBin=hMass->GetBinLowEdge(minBinSig);
 325   Double_t maxMassSigBin=hMass->GetBinLowEdge(maxBinSig)+hMass->GetBinWidth(maxBinSig);
 326   printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 327   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 328   Int_t minBinBkgLow=hMinBin;
 329   Int_t maxBinBkgLow=hMass->FindBin(maxMassBkgLow);
 330   Double_t minMassBkgLowBin=hmin;
 331   Double_t maxMassBkgLowBin=hMass->GetBinLowEdge(maxBinBkgLow)+hMass->GetBinWidth(maxBinBkgLow);
 332   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 333   Int_t minBinBkgHi=hMass->FindBin(minMassBkgHi);
 334   Int_t maxBinBkgHi=hMaxBin;
 335   Double_t minMassBkgHiBin=hMass->GetBinLowEdge(minBinBkgHi);
 336   Double_t maxMassBkgHiBin=hmax;
 337   printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 338   Double_t bkgSig=fB2->Integral(minMassSigBin,maxMassSigBin);
 339   Double_t bkgLow=fB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 340   Double_t bkgHi=fB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 341   printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 342   if(c1){
 343     TBox* bleft=new TBox(minMassBkgLowBin,0.,maxMassBkgLowBin,hMass->GetMaximum());
 344     bleft->SetFillColor(kRed+1);
 345     bleft->SetFillStyle(3002);
 346     bleft->Draw();
 347     TBox* bright=new TBox(minMassBkgHiBin,0.,maxMassBkgHiBin,hMass->GetMaximum());
 348     bright->SetFillColor(kBlue+1);
 349     bright->SetFillStyle(3002);
 350     bright->Draw();
 351     TBox* bsig=new TBox(minMassSigBin,0.,maxMassSigBin,hMass->GetMaximum()*2);
 352     bsig->SetFillColor(1);
 353     bsig->SetFillStyle(3002);
 354     bsig->Draw();
 355   }
 356   TH1F* hCos2PhiBkgLo=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgLoBin%d",iFinalPtBin),minBinBkgLow,maxBinBkgLow);
 357   TH1F* hCos2PhiBkgHi=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgHiBin%d",iFinalPtBin),minBinBkgHi,maxBinBkgHi);
 358   TH1F* hCos2PhiSigReg=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgSigBin%d",iFinalPtBin),minBinSig,maxBinSig);
 359
 360   hCos2PhiBkgLo->Rebin(rebin);
 361   hCos2PhiBkgHi->Rebin(rebin);
 362   hCos2PhiSigReg->Rebin(rebin);
 363   hCos2PhiSigReg->SetLineWidth(2);
 364   hCos2PhiBkgLo->SetLineWidth(2);
 365   hCos2PhiBkgHi->SetLineWidth(2);
 366   hCos2PhiBkgLo->SetLineColor(kRed+1);
 367   hCos2PhiBkgHi->SetLineColor(kBlue+1);
 368
 369   TH1F* hCos2PhiBkgLoScal=(TH1F*)hCos2PhiBkgLo->Clone(Form("hCos2PhiBkgLoScalBin%d",iFinalPtBin));
 370   hCos2PhiBkgLoScal->Scale(bkgSig/bkgLow);
 371   TH1F* hCos2PhiBkgHiScal=(TH1F*)hCos2PhiBkgHi->Clone(Form("hCos2PhiBkgHiScalBin%d",iFinalPtBin));
 372   hCos2PhiBkgHiScal->Scale(bkgSig/bkgHi);
 373   hCos2PhiBkgLoScal->SetLineWidth(2);
 374   hCos2PhiBkgHiScal->SetLineWidth(2);
 375   hCos2PhiBkgLoScal->SetLineColor(kRed+1);
 376   hCos2PhiBkgHiScal->SetLineColor(kBlue+1);
 377   TH1F* hCos2PhiBkgAver=0x0;
 378   if(optBkg==0){
 379     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 380     hCos2PhiBkgAver->Add(hCos2PhiBkgHiScal);
 381     hCos2PhiBkgAver->Scale(0.5);
 382   }else if(optBkg==-1){
 383     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 384   }else{
 385     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgHiScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 386   }
 387   hCos2PhiBkgAver->SetLineWidth(2);
 388   hCos2PhiBkgAver->SetLineColor(kGreen+1);
 389   TH1F* hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone(Form("hCos2PhiSigBin%d",iFinalPtBin));
 390   hCos2PhiSig->Add(hCos2PhiBkgAver,-1.);
 391
 392   TLegend* leg0=new TLegend(0.3,0.6,0.75,0.89);
 393   TPaveText* t0= new TPaveText(0.15,0.70,0.45,0.89,"NDC");
 394   if(c1){
 395     c1->cd(3);
 396     hCos2PhiSigReg->Draw();
 397     hCos2PhiBkgLoScal->Draw("same");
 398     hCos2PhiBkgHiScal->Draw("same");
 399     hCos2PhiBkgAver->Draw("same");
 400     leg0->SetFillColor(0);
 401     TLegendEntry* ent=leg0->AddEntry(hCos2PhiSigReg,"Signal region","L");
 402     ent->SetTextColor(hCos2PhiSigReg->GetLineColor());
 403     ent=leg0->AddEntry(hCos2PhiBkgLoScal,"Left side band","L");
 404     ent->SetTextColor(hCos2PhiBkgLoScal->GetLineColor());
 405     ent=leg0->AddEntry(hCos2PhiBkgHiScal,"Right side band","L");
 406     ent->SetTextColor(hCos2PhiBkgHiScal->GetLineColor());
 407     ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","L");
 408     ent->SetTextColor(hCos2PhiBkgAver->GetLineColor());
 409     leg0->Draw();
 410     c1->cd(4);
 411     hCos2PhiSig->Draw("EP");
 412     t0->SetFillColor(0);
 413     t0->AddText(Form("v2=%.3f+-%.3f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError()));
 414     t0->Draw();
 415
 416     c1->Update();
 417   }
 418   if(!c1){
 419     delete leg0;
 420     delete t0;
 421     delete hCos2PhiBkgLo;
 422     delete hCos2PhiBkgHi;
 423     delete hCos2PhiSigReg;
 424     delete hCos2PhiBkgLoScal;
 425     delete hCos2PhiBkgHiScal;
 426     delete hCos2PhiBkgAver;
 427   }
 428   printf("Signal from mass fitter = %f  Signal from subracted histo= %f\n",
 429          signalFromFit,hCos2PhiSig->Integral());
 430   return hCos2PhiSig;
 431 }
 432
 433
 434 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors, Bool_t useConst){
 435
 436   Int_t npars=fSB->GetNpar();
 437   Double_t sigmaSB=fSB->GetParameter(npars-1);
 438   Double_t massSB=fSB->GetParameter(npars-2);
 439   Double_t integr=fSB->GetParameter(npars-3);
 440   Double_t sOverAll=(fSB->Eval(massSB)-fB2->Eval(massSB))/fSB->Eval(massSB);
 441   printf("mass=%f  S+B=%f   bkg=%f S/(S+B)=%f\n",massSB,fSB->Eval(massSB),fB2->Eval(massSB),sOverAll);
 442   printf("Number of parameters: %d. Signal params: %f %f %f\n",npars,massSB,sigmaSB,integr);
 443   if(optErrors==1){
 444     sOverAll+=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
 445   }
 446   else if(optErrors==-1){
 447     sOverAll-=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
 448   }
 449
 450   Int_t nbinsmass=hMassDphi->GetNbinsY();
 451   Double_t minmass=hMassDphi->GetYaxis()->GetXmin();
 452   Double_t maxmass=hMassDphi->GetYaxis()->GetXmax();
 453
 454   TF1* fSig=new TF1(Form("fSig%d",iFinalPtBin),"[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./[2]/[2])",minmass,maxmass);
 455   fSig->SetParameters(integr,massSB,sigmaSB);
 456
 457   TH1F* hAverCos2Phi=new TH1F(Form("hAverCos2PhiBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 458   TH1F* hFractionSig=new TH1F(Form("hFractionSigBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 459   TH1F* hFractionBkg=new TH1F(Form("hFractionBkgBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 460
 461   for(Int_t iBin=1; iBin<=hMassDphi->GetNbinsY(); iBin++){
 462     TH1F* htemp=(TH1F*)hMassDphi->ProjectionX("htemp",iBin,iBin);
 463     hAverCos2Phi->SetBinContent(iBin,htemp->GetMean());
 464     hAverCos2Phi->SetBinError(iBin,htemp->GetMeanError());
 465     Double_t sig=fSig->Eval(hFractionSig->GetBinCenter(iBin));
 466     Double_t bkg=fB2->Eval(hFractionSig->GetBinCenter(iBin));
 467     if(bkg<1 && sig<1){
 468       hFractionSig->SetBinContent(iBin,0.);
 469       hFractionSig->SetBinError(iBin,0.);
 470       hFractionBkg->SetBinContent(iBin,1.);
 471       hFractionBkg->SetBinError(iBin,0.);
 472     }else{
 473       Double_t fracs=sig/(sig+bkg);
 474       Double_t fracb=bkg/(sig+bkg);
 475       Double_t efracs=0.;//TMath::Sqrt(fracs*(1.-fracs)/(sig+bkg));
 476       Double_t efracb=0.;//TMath::Sqrt(fracb*(1.-fracb)/(sig+bkg));
 477
 478       hFractionSig->SetBinContent(iBin,fracs);
 479       hFractionSig->SetBinError(iBin,efracs);
 480       hFractionBkg->SetBinContent(iBin,fracb);
 481       hFractionBkg->SetBinError(iBin,efracb);
 482     }
 483     delete htemp;
 484   }
 485
 486   TF1* fv2=0x0;
 487   if(useConst){
 488     fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMass,minmass,maxmass,5);
 489   }else{
 490     fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMassLin,minmass,maxmass,6);
 491     fv2->SetParameter(5,0.);
 492   }
 493   fv2->SetParameter(0,sOverAll);
 494   fv2->SetParameter(1,massSB);
 495   fv2->SetParameter(2,sigmaSB);
 496   fv2->SetParameter(3,0.2);
 497   fv2->SetParameter(4,0.2);
 498   fv2->FixParameter(0,sOverAll);
 499   fv2->FixParameter(1,massSB);
 500   fv2->FixParameter(2,sigmaSB);
 501
 502   if((hAverCos2Phi->GetNbinsX()%rebin)==0){
 503     hAverCos2Phi->Rebin(rebin);
 504     hAverCos2Phi->Scale(1./(Double_t)rebin);
 505   }
 506
 507   if(c2){
 508     c2->Divide(2,2);
 509     c2->cd(1);
 510     hMassDphi->Draw("colz");
 511     c2->cd(2);
 512     hMass->Rebin(2);
 513     hMass->SetMinimum(0.);
 514     hMass->SetMarkerStyle(20);
 515     hMass->Draw("E");
 516     fSB->Draw("same");
 517     fSig->Draw("same");
 518     fB2->Draw("same");
 519     c2->cd(3);
 520     hFractionSig->SetMaximum(1.2);
 521     hFractionSig->Draw();
 522     hFractionSig->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 523     hFractionSig->GetYaxis()->SetTitle("Fraction");
 524     hFractionBkg->SetLineColor(2);
 525     hFractionBkg->Draw("same");
 526     TLegend* leg1=new TLegend(0.15,0.15,0.35,0.35);
 527     leg1->SetFillColor(0);
 528     TLegendEntry* ent=leg1->AddEntry(hFractionSig,"S/(S+B)","L");
 529     ent->SetTextColor(hFractionSig->GetLineColor());
 530     ent=leg1->AddEntry(hFractionBkg,"B/(S+B)","L");
 531     ent->SetTextColor(hFractionBkg->GetLineColor());
 532     leg1->Draw();
 533     c2->cd(4);
 534     hAverCos2Phi->Fit(fv2);
 535     fv2->DrawCopy("same");
 536     hAverCos2Phi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 537     hAverCos2Phi->GetYaxis()->SetTitle("v_2^{obs}");
 538     TPaveText* t1= new TPaveText(0.55,0.70,0.89,0.89,"NDC");
 539     t1->SetFillColor(0);
 540     t1->AddText(Form("v2sig=%.3f+-%.3f\n",fv2->GetParameter(3),fv2->GetParError(3)));
 541     if(useConst){
 542       t1->AddText(Form("v2bkg=%.3f+-%.3f\n",fv2->GetParameter(4),fv2->GetParError(4)));
 543     }else{
 544       t1->AddText(Form("v2bkg=(%.3f+-%.3f) + (%.3g+-%.3g)*mass\n",fv2->GetParameter(4),fv2->GetParError(4),fv2->GetParameter(5),fv2->GetParError(5)));
 545     }
 546     t1->Draw();
 547     c2->Update();
 548   }else{
 549     hAverCos2Phi->Fit(fv2);
 550   }
 551   return fv2;
 552 }
 553
 554
 555 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh){
 556   // Event plane resolution syst err (from wide centrality bin
 557   TH1F* hResolSubAB=0x0;
 558   Double_t xmin=1.;
 559   Double_t xmax=-1.;
 560   TGraphAsymmErrors* grSingle=new TGraphAsymmErrors(0);
 561   TGraphAsymmErrors* grInteg=new TGraphAsymmErrors(0);
 562   Int_t iPt=0;
 563   for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){
 564     TString hisnameEP=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+5);
 565     TH1F* hResolSubABsing=(TH1F*)lst->FindObject(hisnameEP.Data());
 566     Double_t resolFull=AliEventPlaneResolution::GetFullEvResol(hResolSubABsing,1);
 567     Double_t resolFullmin=AliEventPlaneResolution::GetFullEvResolLowLim(hResolSubABsing,1);
 568     Double_t resolFullmax=AliEventPlaneResolution::GetFullEvResolHighLim(hResolSubABsing,1);
 569     grSingle->SetPoint(iPt,iHisC+2.5,resolFull);
 570     grSingle->SetPointEXlow(iPt,2.5);
 571     grSingle->SetPointEXhigh(iPt,2.5);
 572     grSingle->SetPointEYlow(iPt,resolFullmax-resolFull);
 573     grSingle->SetPointEYhigh(iPt,resolFull-resolFullmin);
 574     ++iPt;
 575     if(resolFullmin<xmin) xmin=resolFullmin;
 576     if(resolFullmax>xmax) xmax=resolFullmax;
 577     if(iHisC==minCent){
 578       hResolSubAB=(TH1F*)hResolSubABsing->Clone("hResolSubAB");
 579     }else{
 580       hResolSubAB->Add(hResolSubABsing);
 581     }
 582     printf("Adding histogram %s\n",hisnameEP.Data());
 583   }
 584   rcflow=xmin;
 585   rcfhigh=xmax;
 586
 587   Double_t resolSub=AliEventPlaneResolution::GetSubEvResol(hResolSubAB);
 588   printf("Resolution on sub event  = %.4f\n",resolSub);
 589   Double_t chisub=AliEventPlaneResolution::FindChi(resolSub,1);
 590   printf("Chi Subevent             = %.4f\n",chisub);
 591   Double_t chifull=chisub*TMath::Sqrt(2);
 592   printf("Chi Full Event           = %.4f\n",chifull);
 593   Double_t resolFull=AliEventPlaneResolution::GetFullEvResol(hResolSubAB,1);
 594   Double_t resolFullmin=AliEventPlaneResolution::GetFullEvResolLowLim(hResolSubAB,1);
 595   Double_t resolFullmax=AliEventPlaneResolution::GetFullEvResolHighLim(hResolSubAB,1);
 596
 597   AliEventPlaneResolution *resol=new AliEventPlaneResolution(1);
 598   resol->SetSubEventHisto(hResolSubAB);
 599   Double_t resolFull2=resol->GetFullEvResol();
 600   printf("Resolution on full event = %.4f %.4f\n",resolFull,resolFull2);
 601
 602   grInteg->SetPoint(0,40.,resolFull);
 603   grInteg->SetPointEXlow(0,10);
 604   grInteg->SetPointEXhigh(0,10.);
 605   grInteg->SetPointEYlow(0,resolFullmax-resolFull);
 606   grInteg->SetPointEYhigh(0,resolFull-resolFullmin);
 607
 608
 609   TCanvas* cEP=new TCanvas("cEP","EvPlaneResol");
 610   cEP->Divide(1,2);
 611   cEP->cd(1);
 612   hResolSubAB->Draw();
 613   TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
 614   tres->SetNDC();
 615   tres->Draw();
 616   cEP->cd(2);
 617   grSingle->SetMarkerStyle(20);
 618   grInteg->SetMarkerColor(kRed+1);
 619   grInteg->SetLineColor(kRed+1);
 620   grInteg->SetMarkerStyle(29);
 621   grSingle->Draw("AP");
 622   grSingle->GetXaxis()->SetTitle("Centrality Percentile");
 623   grSingle->GetYaxis()->SetTitle("Resolution Correction Factor");
 624   grInteg->Draw("Psame");
 625
 626   return resolFull;
 627
 628 }
 629
 630
 631 void LoadMassHistos(TList* lst, TH2F** hMassDphi){
 632
 633   for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){
 634     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 635       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){
 636         TString hisname=Form("hMc2phi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+5);
 637         TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
 638         if(hMassDphi[iFinalPtBin]==0x0){
 639           hMassDphi[iFinalPtBin]=(TH2F*)htmp->Clone(Form("hMassCos2DphiBin%d",iFinalPtBin));
 640         }else{
 641           hMassDphi[iFinalPtBin]->Add(htmp);
 642         }
 643         printf("Adding histogram %s to PtBin %d\n",hisname.Data(),iFinalPtBin);
 644       }
 645     }
 646   }
 647 }
 648
 649 Bool_t DefinePtBins(TDirectoryFile* df){
 650   AliRDHFCutsD0toKpi *d0cuts=(AliRDHFCutsD0toKpi*)df->Get("D0toKpiCuts");
 651   Int_t nPtBinsCuts=d0cuts->GetNPtBins();
 652   printf("Number of pt bins for cut object = %d\n",nPtBinsCuts);
 653   Float_t *ptlimsCuts=d0cuts->GetPtBinLimits();
 654   for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f \n",iPt,ptlimsCuts[iPt]);
 655   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
 656     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 657       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[iFinalPtBin])<0.0001){
 658         minPtBin[iFinalPtBin]=iPtCuts;
 659         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts;
 660       }
 661     }
 662     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[nFinalPtBins])<0.0001) maxPtBin[nFinalPtBins-1]=iPtCuts;
 663  }
 664
 665   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 666     printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
 667     if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
 668   }
 669   return kTRUE;
 670 }
 671
 672
 673 void SystForSideBands(){
 674   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 675
 676   TFile* filcuts=new TFile(filcutsname.Data());
 677   TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
 678   Bool_t binOK=DefinePtBins(dfcuts);
 679   if(!binOK) return;
 680
 681   TFile* fil=new TFile(filname.Data());
 682   TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
 683   TList* lst=(TList*)df->Get(listname.Data());
 684   if(!lst){
 685     printf("ERROR: list %s not found in file\n",listname.Data());
 686     return;
 687   }
 688
 689   Double_t rcfmin,rcfmax;
 690   Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
 691
 692   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
 693   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 694     hMassDphi[iFinalPtBin]=0x0;
 695   }
 696   LoadMassHistos(lst, hMassDphi);
 697
 698   Int_t nSteps=21;
 699
 700   TGraphErrors** gSystSigRange=new TGraphErrors*[nFinalPtBins];
 701   TGraphErrors** gSystBkgRange=new TGraphErrors*[nFinalPtBins];
 702   TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
 703   TGraphErrors** gSystWhichSide=new TGraphErrors*[nFinalPtBins];
 704
 705   Double_t min1[nFinalPtBins],max1[nFinalPtBins];
 706   Double_t min123[nFinalPtBins],max123[nFinalPtBins];
 707   Double_t min2[nFinalPtBins],max2[nFinalPtBins];
 708   Double_t min3[nFinalPtBins],max3[nFinalPtBins];
 709   Double_t min4[nFinalPtBins],max4[nFinalPtBins];
 710
 711   Double_t sigmaRangeForSigOrig=sigmaRangeForSig;
 712   Double_t sigmaRangeForBkgOrig=sigmaRangeForBkg;
 713
 714
 715   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 716     printf("**************** BIN %d ******************\n",iFinalPtBin);
 717
 718     Int_t rebinHistoSideBandsOrig=rebinHistoSideBands[iFinalPtBin];
 719
 720     gSystSigRange[iFinalPtBin]=new TGraphErrors(0);
 721     gSystSigRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Signal Region Ptbin %d",iFinalPtBin));
 722
 723     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 724
 725     Bool_t out=FitMassSpectrum(hMass,0x0);
 726     if(!out) continue;
 727
 728     min1[iFinalPtBin]=99999.;
 729     max1[iFinalPtBin]=-99999.;
 730     min123[iFinalPtBin]=99999.;
 731     max123[iFinalPtBin]=-99999.;
 732     for(Int_t iStep=0; iStep<nSteps; iStep++){
 733       Int_t index=iFinalPtBin*nSteps+iStep;
 734       sigmaRangeForSig=1.5+0.1*iStep;
 735       sigmaRangeForBkg=sigmaRangeForBkgOrig;
 736       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 737       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 738       Double_t v2obsM1=hCos2PhiSig->GetMean();
 739       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 740       delete hCos2PhiSig;
 741       Double_t v2M1=v2obsM1/resolFull;
 742       Double_t errv2M1=errv2obsM1/resolFull;
 743       gSystSigRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForSig,v2M1);
 744       gSystSigRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
 745       if(v2M1>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M1;
 746       if(v2M1<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M1;
 747       if(sigmaRangeForSig>=2. && sigmaRangeForSig<=3){
 748         if(v2M1>max123[iFinalPtBin]) max123[iFinalPtBin]=v2M1;
 749         if(v2M1<min123[iFinalPtBin]) min123[iFinalPtBin]=v2M1;
 750       }
 751     }
 752
 753     min2[iFinalPtBin]=99999.;
 754     max2[iFinalPtBin]=-99999.;
 755     gSystBkgRange[iFinalPtBin]=new TGraphErrors(0);
 756     gSystBkgRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Bkg Region Ptbin %d",iFinalPtBin));
 757     for(Int_t iStep=0; iStep<nSteps; iStep++){
 758       Int_t index=nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
 759       sigmaRangeForSig=sigmaRangeForSigOrig;
 760       sigmaRangeForBkg=4.+0.1*iStep;
 761       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 762       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 763       Double_t v2obsM1=hCos2PhiSig->GetMean();
 764       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 765       delete hCos2PhiSig;
 766       Double_t v2M1=v2obsM1/resolFull;
 767       Double_t errv2M1=errv2obsM1/resolFull;
 768       gSystBkgRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForBkg,v2M1);
 769       gSystBkgRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
 770       if(v2M1>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M1;
 771       if(v2M1<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M1;
 772     }
 773
 774     min3[iFinalPtBin]=99999.;
 775     max3[iFinalPtBin]=-99999.;
 776     gSystRebin[iFinalPtBin]=new TGraphErrors(0);
 777     gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin Ptbin %d",iFinalPtBin));
 778     Int_t nPts=0;
 779     for(Int_t iStep=0; iStep<nSteps; iStep++){
 780       Int_t index=2*nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
 781       sigmaRangeForSig=sigmaRangeForSigOrig;
 782       sigmaRangeForBkg=sigmaRangeForBkgOrig;
 783       rebinHistoSideBands[iFinalPtBin]=iStep+1;
 784       if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistoSideBands[iFinalPtBin]!=0) continue;
 785       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 786       Double_t v2obsM1=hCos2PhiSig->GetMean();
 787       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 788       delete hCos2PhiSig;
 789       Double_t v2M1=v2obsM1/resolFull;
 790       Double_t errv2M1=errv2obsM1/resolFull;
 791       gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistoSideBands[iFinalPtBin],v2M1);
 792       gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
 793       if(v2M1>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M1;
 794       if(v2M1<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M1;
 795       ++nPts;
 796     }
 797
 798     min4[iFinalPtBin]=99999.;
 799     max4[iFinalPtBin]=-99999.;
 800     gSystWhichSide[iFinalPtBin]=new TGraphErrors(0);
 801     gSystWhichSide[iFinalPtBin]->SetTitle(Form("v2 vs. WhichSide Ptbin %d",iFinalPtBin));
 802     nPts=0;
 803     for(Int_t iStep=-1; iStep<=1; iStep++){
 804       Int_t index=3*nSteps*nFinalPtBins+2+iStep;
 805       sigmaRangeForSig=sigmaRangeForSigOrig;
 806       sigmaRangeForBkg=sigmaRangeForBkgOrig;
 807       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 808       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0,iStep);
 809       Double_t v2obsM1=hCos2PhiSig->GetMean();
 810       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 811       delete hCos2PhiSig;
 812       Double_t v2M1=v2obsM1/resolFull;
 813       Double_t errv2M1=errv2obsM1/resolFull;
 814       gSystWhichSide[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M1);
 815       gSystWhichSide[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
 816       if(v2M1>max4[iFinalPtBin]) max4[iFinalPtBin]=v2M1;
 817       if(v2M1<min4[iFinalPtBin]) min4[iFinalPtBin]=v2M1;
 818       ++nPts;
 819     }
 820     rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 821   }
 822
 823   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 824     printf("------ Pt Bin %d ------\n",iFinalPtBin);
 825     printf("Range of values for sig variation = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
 826     printf("           (limited to 2-3 sigma) = %f %f\n",min123[iFinalPtBin],max123[iFinalPtBin]);
 827     printf("Range of values for bkg variation = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
 828     printf("Range of values for rebin = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
 829     printf("Range of values for whichside = %f %f\n",min4[iFinalPtBin],max4[iFinalPtBin]);
 830     Float_t minenv=TMath::Min(min123[iFinalPtBin],TMath::Min(min2[iFinalPtBin],min4[iFinalPtBin]));
 831     Float_t maxenv=TMath::Max(max123[iFinalPtBin],TMath::Max(max2[iFinalPtBin],max4[iFinalPtBin]));
 832     printf(" ---> Envelope                = %f %f\n",minenv,maxenv);
 833   }
 834
 835   TCanvas* cs1=new TCanvas("cs1");
 836   cs1->Divide(2,2);
 837   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 838     cs1->cd(iFinalPtBin+1);
 839     gSystSigRange[iFinalPtBin]->SetMarkerStyle(20);
 840     gSystSigRange[iFinalPtBin]->Draw("AP");
 841     gSystSigRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaSignal");
 842     gSystSigRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
 843   }
 844
 845   TCanvas* cs2=new TCanvas("cs2");
 846   cs2->Divide(2,2);
 847   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 848     cs2->cd(iFinalPtBin+1);
 849     gSystBkgRange[iFinalPtBin]->SetMarkerStyle(20);
 850     gSystBkgRange[iFinalPtBin]->Draw("AP");
 851     gSystBkgRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaBackground");
 852     gSystBkgRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
 853   }
 854
 855   TCanvas* cs3=new TCanvas("cs3");
 856   cs3->Divide(2,2);
 857   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 858     cs3->cd(iFinalPtBin+1);
 859     gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
 860     gSystRebin[iFinalPtBin]->Draw("AP");
 861     gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
 862     gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
 863   }
 864
 865   TCanvas* cs4=new TCanvas("cs4");
 866   cs4->Divide(2,2);
 867   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 868     cs4->cd(iFinalPtBin+1);
 869     gSystWhichSide[iFinalPtBin]->SetMarkerStyle(20);
 870     gSystWhichSide[iFinalPtBin]->Draw("AP");
 871     gSystWhichSide[iFinalPtBin]->GetXaxis()->SetTitle("Side band used");
 872     gSystWhichSide[iFinalPtBin]->GetYaxis()->SetTitle("v2");
 873   }
 874 }
 875
 876 void SystForFitv2Mass(){
 877   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 878
 879   TFile* filcuts=new TFile(filcutsname.Data());
 880   TDirectoryFile* dfcuts=(TDirectoryFile*)filcuts->Get("PWG3_D2H_HFv2");
 881   Bool_t binOK=DefinePtBins(dfcuts);
 882   if(!binOK) return;
 883
 884   TFile* fil=new TFile(filname.Data());
 885   TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
 886   TList* lst=(TList*)df->Get(listname.Data());
 887   if(!lst){
 888     printf("ERROR: list %s not found in file\n",listname.Data());
 889     return;
 890   }
 891
 892   Double_t rcfmin,rcfmax;
 893   Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
 894
 895   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
 896   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 897     hMassDphi[iFinalPtBin]=0x0;
 898   }
 899   LoadMassHistos(lst, hMassDphi);
 900
 901   Int_t nSteps=11;
 902
 903   TGraphErrors** gSystParamErr=new TGraphErrors*[nFinalPtBins];
 904   TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
 905   TGraphErrors** gSystLinConst=new TGraphErrors*[nFinalPtBins];
 906
 907   Double_t min1[nFinalPtBins],max1[nFinalPtBins];
 908   Double_t min2[nFinalPtBins],max2[nFinalPtBins];
 909   Double_t min3[nFinalPtBins],max3[nFinalPtBins];
 910
 911
 912   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 913     printf("**************** BIN %d ******************\n",iFinalPtBin);
 914
 915     Int_t rebinHistov2MassOrig=rebinHistov2Mass[iFinalPtBin];
 916     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 917
 918     Bool_t out=FitMassSpectrum(hMass,0x0);
 919     if(!out) continue;
 920
 921     min1[iFinalPtBin]=99999.;
 922     max1[iFinalPtBin]=-99999.;
 923     gSystRebin[iFinalPtBin]=new TGraphErrors(0);
 924     gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin PtBin %d",iFinalPtBin));
 925     Int_t nPts=0;
 926     for(Int_t iStep=0; iStep<nSteps; iStep++){
 927       Int_t index=iStep;
 928       rebinHistov2Mass[iFinalPtBin]=iStep+1;
 929       if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistov2Mass[iFinalPtBin]!=0) continue;
 930       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,useConstantvsBkgVsMass);
 931       Double_t v2obsM2=fv2->GetParameter(3);
 932       Double_t errv2obsM2=fv2->GetParError(3);
 933       delete fv2;
 934       Double_t v2M2=v2obsM2/resolFull;
 935       Double_t errv2M2=errv2obsM2/resolFull;
 936       gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistov2Mass[iFinalPtBin],v2M2);
 937       gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
 938       if(v2M2>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M2;
 939       if(v2M2<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M2;
 940       ++nPts;
 941     }
 942     rebinHistov2Mass[iFinalPtBin]=rebinHistov2MassOrig;
 943
 944     min2[iFinalPtBin]=99999.;
 945     max2[iFinalPtBin]=-99999.;
 946     gSystParamErr[iFinalPtBin]=new TGraphErrors(0);
 947     gSystParamErr[iFinalPtBin]->SetTitle(Form("v2 vs. ParamErr PtBin %d",iFinalPtBin));
 948     nPts=0;
 949     for(Int_t iStep=-1; iStep<=1; iStep++){
 950       Int_t index=nSteps*2+iStep;
 951       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,iStep,useConstantvsBkgVsMass);
 952       Double_t v2obsM2=fv2->GetParameter(3);
 953       Double_t errv2obsM2=fv2->GetParError(3);
 954       delete fv2;
 955       Double_t v2M2=v2obsM2/resolFull;
 956       Double_t errv2M2=errv2obsM2/resolFull;
 957       gSystParamErr[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M2);
 958       gSystParamErr[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
 959       if(v2M2>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M2;
 960       if(v2M2<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M2;
 961       ++nPts;
 962     }
 963
 964     min3[iFinalPtBin]=99999.;
 965     max3[iFinalPtBin]=-99999.;
 966     gSystLinConst[iFinalPtBin]=new TGraphErrors(0);
 967     gSystLinConst[iFinalPtBin]->SetTitle(Form("v2 LinVsConst Bin%d",iFinalPtBin));
 968     nPts=0;
 969     for(Int_t iStep=0; iStep<=1; iStep++){
 970       Int_t index=nSteps*3+iStep;
 971       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,iStep);
 972       Double_t v2obsM2=fv2->GetParameter(3);
 973       Double_t errv2obsM2=fv2->GetParError(3);
 974       delete fv2;
 975       Double_t v2M2=v2obsM2/resolFull;
 976       Double_t errv2M2=errv2obsM2/resolFull;
 977       gSystLinConst[iFinalPtBin]->SetPoint(nPts,1.-(Double_t)iStep,v2M2);
 978       gSystLinConst[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
 979       if(v2M2>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M2;
 980       if(v2M2<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M2;
 981       ++nPts;
 982     }
 983
 984
 985   }
 986
 987   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 988     printf("------ Pt Bin %d ------\n",iFinalPtBin);
 989     printf("Range of values for rebin = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
 990     printf("Range of values for par err = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
 991     printf("Range of values for lin const = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
 992     Float_t minenv=TMath::Min(min2[iFinalPtBin],min3[iFinalPtBin]);
 993     Float_t maxenv=TMath::Max(max2[iFinalPtBin],max3[iFinalPtBin]);
 994     printf(" ---> Envelope                = %f %f\n",minenv,maxenv);
 995   }
 996
 997   TCanvas* cs1=new TCanvas("cs1");
 998   cs1->Divide(2,2);
 999   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1000     cs1->cd(iFinalPtBin+1);
1001     gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
1002     gSystRebin[iFinalPtBin]->Draw("AP");
1003     gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
1004     gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1005   }
1006
1007   TCanvas* cs2=new TCanvas("cs2");
1008   cs2->Divide(2,2);
1009   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1010     cs2->cd(iFinalPtBin+1);
1011     gSystParamErr[iFinalPtBin]->SetMarkerStyle(20);
1012     gSystParamErr[iFinalPtBin]->Draw("AP");
1013     gSystParamErr[iFinalPtBin]->GetXaxis()->SetTitle("Error on Signal yield");
1014     gSystParamErr[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1015   }
1016
1017   TCanvas* cs3=new TCanvas("cs3");
1018   cs3->Divide(2,2);
1019   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1020     cs3->cd(iFinalPtBin+1);
1021     gSystLinConst[iFinalPtBin]->SetMarkerStyle(20);
1022     gSystLinConst[iFinalPtBin]->Draw("AP");
1023     gSystLinConst[iFinalPtBin]->GetXaxis()->SetLimits(-0.1,1.1);
1024     gSystLinConst[iFinalPtBin]->GetXaxis()->SetTitle("Const/Linear v2 of background");
1025     gSystLinConst[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1026   }
1027 }