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