PWGHF/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 #include "AliAODRecoDecayHF.h"
  23 #include "AliRDHFCuts.h"
  24 #include "AliRDHFCutsDplustoKpipi.h"
  25 #include "AliRDHFCutsDStartoKpipi.h"
  26 #include "AliRDHFCutsD0toKpi.h"
  27 #include "AliHFMassFitter.h"
  28 #include "AliEventPlaneResolutionHandler.h"
  29 #include "AliVertexingHFUtils.h"
  30 #endif
  31
  32
  33 // Common variables: to be configured by the user
  34 TString filename="AnalysisResults_train63.root";
  35 TString suffix="v2Dplus3050Cut4upcutPIDTPC";
  36 TString partname="Dplus";
  37 Int_t minCent=30;
  38 Int_t maxCent=50;
  39 //evPlane flag from AliEventPlaneResolutionHandler:
  40 //kTPCFullEta, kTPCPosEta,kVZERO,kVZEROA,kVZEROC
  41 Int_t evPlane=AliEventPlaneResolutionHandler::kTPCPosEta;
  42 //resolution flag fromAliEventPlaneResolutionHandler:
  43 //kTwoRandSub,kTwoChargeSub,kTwoEtaSub,kThreeSub,kThreeSubTPCGap
  44 Int_t evPlaneRes=AliEventPlaneResolutionHandler::kTwoEtaSub;
  45 Bool_t useNcollWeight=kFALSE;
  46
  47 const Int_t nFinalPtBins=4;
  48 Double_t ptlims[nFinalPtBins+1]={3.,4.,6.,8.,12.};
  49 Double_t sigmaRangeForSig=2.5;
  50 Double_t sigmaRangeForBkg=4.5;
  51 Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2,2};
  52 Bool_t useConstantvsBkgVsMass=kFALSE;
  53 Int_t rebinHistov2Mass[nFinalPtBins]={2,2,2,2};
  54 Int_t factor4refl=0;
  55 Int_t minPtBin[nFinalPtBins]={-1,-1,-1,-1};
  56 Int_t maxPtBin[nFinalPtBins]={-1,-1,-1,-1};
  57 Bool_t saveAllCanvas=kFALSE;
  58
  59 // systematic errors for D+, 30-50 TPC eventplane
  60 Double_t systErrMeth1[nFinalPtBins]={
  61   (0.247359-0.152503)/2.,
  62   (0.337073-0.290978)/2.,
  63   (0.286888-0.222045)/2.,
  64   (-0.078888+0.178259)/2.
  65 };
  66
  67 Double_t systErrMeth2[nFinalPtBins]={
  68   (0.151019-0.066624)/2.,
  69   (0.272252-0.232044)/2.,
  70   (0.246019-0.143002)/2.,
  71   (-0.026077+0.292956)/2.
  72 };
  73
  74 // systematic errors for 2-5, 5-8 and 8-12 (no run-by-run weights)
  75 /*
  76 Double_t systErrMeth1[nFinalPtBins]={
  77   (0.308-0.169)/2.,
  78   (0.14-0.1)/2.,
  79   (0.04+0.02)/2.
  80 };
  81 Double_t systErrMeth2[nFinalPtBins]={
  82   (0.305-0.252)/2.,
  83   (0.129-0.020)/2.,
  84   (0.101+0.06)/2.
  85 };
  86 */
  87
  88 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights)
  89 /*
  90 Double_t systErrMeth1[nFinalPtBins]={
  91   (0.23-0.10)/2.,
  92   (0.152-0.078)/2.,
  93   (0.161-0.097)/2.
  94 };
  95 Double_t systErrMeth2[nFinalPtBins]={
  96   (0.164-0.097)/2.,
  97   (0.110-0.012)/2.,
  98   (0.131-0.036)/2.
  99 ;
 100 */
 101
 102
 103
 104 // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights, RAA cuts)
 105 /*
 106 Double_t systErrMeth1[nFinalPtBins]={
 107   (0.265-0.122)/2.,
 108   (0.165-0.117)/2.,
 109   (0.238-0.169)/2.
 110 };
 111 Double_t systErrMeth2[nFinalPtBins]={
 112   (0.174-0.135)/2.,
 113   (0.18-0.11)/2.,
 114   (0.311-0.28)/2.
 115 };
 116 */
 117
 118 // output of mass fitter
 119 Int_t hMinBin;
 120 Int_t hMaxBin;
 121 Double_t signalFromFit,esignalFromFit;
 122 Double_t massFromFit,sigmaFromFit;
 123 TF1* fB1=0x0;
 124 TF1* fB2=0x0;
 125 TF1* fSB=0x0;
 126
 127 void LoadMassHistos(TList* lst, TH2F** hMassDphi);
 128 Bool_t DefinePtBins(TDirectoryFile* df);
 129 Double_t v2vsMass(Double_t *x, Double_t *par);
 130 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh);
 131 Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);
 132 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad);
 133 TH1F* DoSideBands(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c1, Int_t optBkg=0);
 134 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors=0, Bool_t useConst=kTRUE);
 135
 136
 137 //______________________________________________________________________________
 138 void Extractv2from2Dhistos(){
 139   // main function: computes v2 with side band and v2(M) methods
 140
 141   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 142
 143   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
 144   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 145
 146   TFile *f = TFile::Open(filename.Data());
 147   if(!f){
 148     printf("file %s not found, please check file name\n",filename.Data());
 149     return;
 150   }
 151   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 152   if(!dir){
 153     printf("Directory %s not found, please check dir name\n",dirname.Data());
 154     return;
 155   }
 156   Bool_t binOK=DefinePtBins(dir);
 157   if(!binOK){
 158     printf("ERROR: mismatch in pt binning\n");
 159     return;
 160   }
 161   TList *lst =(TList*)dir->Get(listname.Data());
 162   if(!lst){
 163     printf("list %s not found in file, please check list name\n",listname.Data());
 164     return;
 165   }
 166
 167   // Double_t rcfmin,rcfmax;
 168   // Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
 169   // Double_t resolSyst=(rcfmax-rcfmin)/2./resolFull;
 170   // printf("Relative Systematic Error on RCF=%f\n",resolSyst);
 171   AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
 172   epres->SetEventPlane(evPlane);
 173   epres->SetResolutionOption(evPlaneRes);
 174   epres->SetUseNcollWeights(useNcollWeight);
 175   Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
 176   Double_t rcfmin=epres->GetEventPlaneResolution(minCent,minCent+2.5);
 177   Double_t rcfmax=epres->GetEventPlaneResolution(maxCent-2.5,maxCent);
 178   Double_t resolSyst=TMath::Abs(rcfmax-rcfmin)/2./resolFull;
 179   delete epres;
 180   printf("Event plane resolution %f\n",resolFull);
 181
 182   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
 183   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 184     hMassDphi[iFinalPtBin]=0x0;
 185   }
 186   LoadMassHistos(lst, hMassDphi);
 187
 188   TCanvas** c1=new TCanvas*[nFinalPtBins];
 189   TCanvas** c2=new TCanvas*[nFinalPtBins];
 190
 191   Double_t v2M1[nFinalPtBins],errv2M1[nFinalPtBins];
 192   Double_t v2M2[nFinalPtBins],errv2M2[nFinalPtBins];
 193
 194   gStyle->SetPalette(1);
 195   gStyle->SetOptTitle(0);
 196
 197   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 198     printf("**************** BIN %d ******************\n",iFinalPtBin);
 199     printf("\n--------- Method 1: Side Bands ----------\n");
 200     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 201
 202     c1[iFinalPtBin]=new TCanvas(Form("cMeth1PtBin%d",iFinalPtBin),Form("Meth1PtBin%d",iFinalPtBin));
 203     c1[iFinalPtBin]->Divide(2,2);
 204     c1[iFinalPtBin]->cd(1);
 205     hMassDphi[iFinalPtBin]->Draw("colz");
 206     c1[iFinalPtBin]->cd(2);
 207     Bool_t out=FitMassSpectrum(hMass,(TPad*)gPad);
 208     if(!out) continue;
 209
 210     TH1F* hCos2PhiSig=DoSideBands(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],c1[iFinalPtBin]);
 211     Double_t v2obsM1=hCos2PhiSig->GetMean();
 212     Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 213     printf("v2obs = %f +- %f\n",v2obsM1,errv2obsM1);
 214     v2M1[iFinalPtBin]=v2obsM1/resolFull;
 215     errv2M1[iFinalPtBin]=errv2obsM1/resolFull;
 216     printf("v2 = %f +- %f\n",v2M1[iFinalPtBin],errv2M1[iFinalPtBin]);
 217
 218     printf("\n--------- Method 2: S/S+B ----------\n");
 219     c2[iFinalPtBin]=new TCanvas(Form("cMeth2PtBin%d",iFinalPtBin),Form("Meth2PtBin%d",iFinalPtBin));
 220     TF1* fv2=DoFitv2vsMass(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],c2[iFinalPtBin],0,useConstantvsBkgVsMass);
 221
 222     Double_t v2obsM2=fv2->GetParameter(3);
 223     Double_t errv2obsM2=fv2->GetParError(3);
 224     printf("v2obs = %f +- %f\n",v2obsM2,errv2obsM2);
 225     v2M2[iFinalPtBin]=v2obsM2/resolFull;
 226     errv2M2[iFinalPtBin]=errv2obsM2/resolFull;
 227     printf("v2 = %f +- %f\n",v2M2[iFinalPtBin],errv2M2[iFinalPtBin]);
 228     if(saveAllCanvas){
 229       c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.root",iFinalPtBin));
 230       c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.root",iFinalPtBin));
 231       c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.eps",iFinalPtBin));
 232       c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.eps",iFinalPtBin));
 233    }
 234   }
 235
 236   printf("\n--------- Summary ------------\n");
 237
 238   TH1F* hv2m1=new TH1F("hv2m1","Side Band subtraction",nFinalPtBins,ptlims);
 239   TH1F* hv2m2=new TH1F("hv2m2","Fit to v2 vs. mass",nFinalPtBins,ptlims);
 240    for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 241     printf("PtBin %d   v2method1 = %f +- %f   v2method2 = %f +-%f\n",iFinalPtBin,
 242            v2M1[iFinalPtBin],errv2M1[iFinalPtBin],
 243            v2M2[iFinalPtBin],errv2M2[iFinalPtBin]
 244            );
 245     hv2m1->SetBinContent(iFinalPtBin+1,v2M1[iFinalPtBin]);
 246     hv2m1->SetBinError(iFinalPtBin+1,errv2M1[iFinalPtBin]);
 247     hv2m2->SetBinContent(iFinalPtBin+1,v2M2[iFinalPtBin]);
 248     hv2m2->SetBinError(iFinalPtBin+1,errv2M2[iFinalPtBin]);
 249   }
 250
 251    TH1F* hSystErr1=(TH1F*)hv2m1->Clone("hSystErr1");
 252    TH1F* hSystErr2=(TH1F*)hv2m2->Clone("hSystErr2");
 253    for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 254      Double_t systRes1=resolSyst*v2M1[iFinalPtBin];
 255      Double_t systRes2=resolSyst*v2M2[iFinalPtBin];
 256      Double_t syste1=TMath::Sqrt(systErrMeth1[iFinalPtBin]*systErrMeth1[iFinalPtBin]+systRes1*systRes1);
 257      Double_t syste2=TMath::Sqrt(systErrMeth2[iFinalPtBin]*systErrMeth2[iFinalPtBin]+systRes2*systRes2);
 258      hSystErr1->SetBinError(iFinalPtBin+1,syste1);
 259      hSystErr2->SetBinError(iFinalPtBin+1,syste2);
 260    }
 261
 262
 263    Double_t maxy=TMath::Max(hv2m2->GetMaximum(),hv2m1->GetMaximum())+0.1;
 264    Double_t miny=TMath::Min(hv2m2->GetMinimum(),hv2m1->GetMinimum())-0.1;
 265    TH2F* hempty=new TH2F("hempty","",10,0.,hv2m1->GetXaxis()->GetXmax()+2.,10,miny,maxy);
 266    hempty->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 267    hempty->GetYaxis()->SetTitle("v_{2}");
 268
 269    TCanvas* cv2=new TCanvas("cv2","v2");
 270    hempty->Draw();
 271    hv2m1->SetMarkerStyle(26);
 272    hSystErr1->SetFillColor(kGray);
 273    hSystErr1->SetFillStyle(3002);
 274    hSystErr1->Draw("e2same");
 275
 276    hv2m2->SetLineColor(kRed+1);
 277    hv2m2->SetMarkerColor(kRed+1);
 278    hv2m2->SetMarkerStyle(20);
 279    hSystErr2->SetFillColor(kRed-9);
 280    hSystErr2->SetFillStyle(3005);
 281    hSystErr2->Draw("e2same");
 282
 283    hv2m1->Draw("same");
 284    hv2m2->Draw("same");
 285
 286    TLegend* leg2=new TLegend(0.5,0.7,0.89,0.89);
 287    leg2->SetFillStyle(0);
 288    TLegendEntry* ent=leg2->AddEntry(hv2m1,"Side Band subtraction","P");
 289    ent->SetTextColor(hv2m1->GetMarkerColor());
 290    ent=leg2->AddEntry(hv2m2,"Fit to v2 vs. mass","P");
 291    ent->SetTextColor(hv2m2->GetMarkerColor());
 292    leg2->Draw();
 293    cv2->Update();
 294    cv2->SaveAs(Form("v2%s-2Dmethods-%s.png",partname.Data(),suffix.Data()));
 295    cv2->SaveAs(Form("v2%s-2Dmethods-%s.eps",partname.Data(),suffix.Data()));
 296
 297    TFile* outfil=new TFile(Form("v2Output2DMeth_%s_%d_%d_%s.root",partname.Data(),minCent,maxCent,suffix.Data()),"recreate");
 298    outfil->cd();
 299    hv2m1->Write();
 300    hv2m2->Write();
 301    hSystErr1->Write();
 302    hSystErr2->Write();
 303    outfil->Close();
 304 }
 305
 306
 307 //______________________________________________________________________________
 308 Double_t v2vsMass(Double_t *x, Double_t *par){
 309   // Fit function for signal+background
 310   // par[0] = S/B at the mass peak
 311   // par[1] = mass
 312   // par[2] = sigma
 313   // par[3] = v2sig
 314   // par[4] = v2back
 315
 316   Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
 317   Double_t fracbkg=1-fracsig;
 318   return par[3]*fracsig+par[4]*fracbkg;
 319 }
 320
 321 //______________________________________________________________________________
 322 Double_t v2vsMassLin(Double_t *x, Double_t *par){
 323   // Fit function for signal+background
 324   // par[0] = S/B at the mass peak
 325   // par[1] = mass
 326   // par[2] = sigma
 327   // par[3] = v2sig
 328   // par[4] = v2back constant
 329   // par[5] = v2back slope
 330
 331   Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
 332   Double_t fracbkg=1-fracsig;
 333   return par[3]*fracsig+(par[4]+par[5]*x[0])*fracbkg;
 334 }
 335
 336 //______________________________________________________________________________
 337 Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
 338
 339   Int_t nMassBins=hMass->GetNbinsX();
 340   hMinBin=3;
 341   hMaxBin=nMassBins-2;
 342   Double_t hmin=hMass->GetBinLowEdge(hMinBin);
 343   Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
 344   Double_t massD;
 345   if(partname.Contains("Dzero")) {
 346     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 347   }else if(partname.Contains("Dplus")){
 348     massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 349   }else if(partname.Contains("Dstar")) {
 350     massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass());
 351   }else{
 352     printf("Wrong particle name\n");
 353     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 354   }
 355
 356   AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
 357   fitter->SetReflectionSigmaFactor(factor4refl);
 358   fitter->SetInitialGaussianMean(massD);
 359   Bool_t out=fitter->MassFitter(0);
 360   if(!out) return kFALSE;
 361   fitter->Signal(sigmaRangeForSig, signalFromFit,esignalFromFit);
 362   massFromFit=fitter->GetMean();
 363   sigmaFromFit=fitter->GetSigma();
 364   fB1=fitter->GetBackgroundFullRangeFunc();
 365   fB2=fitter->GetBackgroundRecalcFunc();
 366   fSB=fitter->GetMassFunc();
 367   if(!fB1) return kFALSE;
 368   if(!fB2) return kFALSE;
 369   if(!fSB) return kFALSE;
 370   if(pad){
 371     fitter->DrawHere(gPad,3.,0.);
 372   }
 373   return kTRUE;
 374 }
 375
 376 //______________________________________________________________________________
 377 TH1F* DoSideBands(Int_t iFinalPtBin,
 378                   TH2F* hMassDphi,
 379                   TH1F* hMass,
 380                   Int_t rebin,
 381                   TCanvas* c1,
 382                   Int_t optBkg){
 383
 384   // Build histo with cos(2*deltaphi) distribution for signal
 385
 386   Double_t hmin=hMass->GetBinLowEdge(hMinBin);
 387   Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
 388
 389   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
 390   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
 391   Int_t minBinSig=hMass->FindBin(minMassSig);
 392   Int_t maxBinSig=hMass->FindBin(maxMassSig);
 393   Double_t minMassSigBin=hMass->GetBinLowEdge(minBinSig);
 394   Double_t maxMassSigBin=hMass->GetBinLowEdge(maxBinSig)+hMass->GetBinWidth(maxBinSig);
 395   printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 396   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
 397   Int_t minBinBkgLow=hMinBin;
 398   Int_t maxBinBkgLow=hMass->FindBin(maxMassBkgLow);
 399   Double_t minMassBkgLowBin=hmin;
 400   Double_t maxMassBkgLowBin=hMass->GetBinLowEdge(maxBinBkgLow)+hMass->GetBinWidth(maxBinBkgLow);
 401   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
 402   Int_t minBinBkgHi=hMass->FindBin(minMassBkgHi);
 403   Int_t maxBinBkgHi=hMaxBin;
 404   Double_t minMassBkgHiBin=hMass->GetBinLowEdge(minBinBkgHi);
 405   Double_t maxMassBkgHiBin=hmax;
 406   printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 407   Double_t bkgSig=fB2->Integral(minMassSigBin,maxMassSigBin);
 408   Double_t bkgLow=fB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 409   Double_t bkgHi=fB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 410   printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 411   if(c1){
 412     TBox* bleft=new TBox(minMassBkgLowBin,0.,maxMassBkgLowBin,hMass->GetMaximum());
 413     bleft->SetFillColor(kRed+1);
 414     bleft->SetFillStyle(3002);
 415     bleft->Draw();
 416     TBox* bright=new TBox(minMassBkgHiBin,0.,maxMassBkgHiBin,hMass->GetMaximum());
 417     bright->SetFillColor(kBlue+1);
 418     bright->SetFillStyle(3002);
 419     bright->Draw();
 420     TBox* bsig=new TBox(minMassSigBin,0.,maxMassSigBin,hMass->GetMaximum()*2);
 421     bsig->SetFillColor(1);
 422     bsig->SetFillStyle(3002);
 423     bsig->Draw();
 424   }
 425   TH1F* hCos2PhiBkgLo=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgLoBin%d",iFinalPtBin),minBinBkgLow,maxBinBkgLow);
 426   TH1F* hCos2PhiBkgHi=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgHiBin%d",iFinalPtBin),minBinBkgHi,maxBinBkgHi);
 427   TH1F* hCos2PhiSigReg=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgSigBin%d",iFinalPtBin),minBinSig,maxBinSig);
 428
 429   printf("Before Rebin11: BkgLo = %f+-%f BkgHi=%f+-%f  SigReg=%f+-%f\n",hCos2PhiBkgLo->GetBinContent(11),hCos2PhiBkgLo->GetBinError(11),
 430          hCos2PhiBkgHi->GetBinContent(11),hCos2PhiBkgHi->GetBinError(11),
 431          hCos2PhiSigReg->GetBinContent(11),hCos2PhiSigReg->GetBinError(11));
 432   printf("Before Rebin12: BkgLo = %f+-%f BkgHi=%f+-%f  SigReg=%f+-%f\n",hCos2PhiBkgLo->GetBinContent(12),hCos2PhiBkgLo->GetBinError(12),
 433          hCos2PhiBkgHi->GetBinContent(12),hCos2PhiBkgHi->GetBinError(12),
 434          hCos2PhiSigReg->GetBinContent(12),hCos2PhiSigReg->GetBinError(12));
 435   Double_t binc=hCos2PhiBkgLo->GetBinCenter(11);
 436
 437   hCos2PhiBkgLo->Rebin(rebin);
 438   hCos2PhiBkgHi->Rebin(rebin);
 439   hCos2PhiSigReg->Rebin(rebin);
 440   hCos2PhiSigReg->SetLineWidth(2);
 441   hCos2PhiBkgLo->SetLineWidth(2);
 442   hCos2PhiBkgHi->SetLineWidth(2);
 443   hCos2PhiBkgLo->SetLineColor(kRed+1);
 444   hCos2PhiBkgHi->SetLineColor(kBlue+1);
 445
 446   Int_t theBinR=hCos2PhiBkgLo->FindBin(binc);
 447   printf("After Rebin %d: BkgLo = %f+-%f BkgHi=%f+-%f  SigReg=%f+-%f\n",theBinR,
 448          hCos2PhiBkgLo->GetBinContent(theBinR),hCos2PhiBkgLo->GetBinError(theBinR),
 449          hCos2PhiBkgHi->GetBinContent(theBinR),hCos2PhiBkgHi->GetBinError(theBinR),
 450          hCos2PhiSigReg->GetBinContent(theBinR),hCos2PhiSigReg->GetBinError(theBinR));
 451
 452   printf("Background Scaling: Lo->Sig= %f Hi->Sig=%f\n",bkgSig/bkgLow,bkgSig/bkgHi);
 453
 454   hCos2PhiBkgLo->Sumw2();
 455   hCos2PhiBkgHi->Sumw2();
 456   TH1F* hCos2PhiBkgLoScal=(TH1F*)hCos2PhiBkgLo->Clone(Form("hCos2PhiBkgLoScalBin%d",iFinalPtBin));
 457   hCos2PhiBkgLoScal->Scale(bkgSig/bkgLow);
 458   TH1F* hCos2PhiBkgHiScal=(TH1F*)hCos2PhiBkgHi->Clone(Form("hCos2PhiBkgHiScalBin%d",iFinalPtBin));
 459   hCos2PhiBkgHiScal->Scale(bkgSig/bkgHi);
 460   hCos2PhiBkgLoScal->SetLineWidth(2);
 461   hCos2PhiBkgHiScal->SetLineWidth(2);
 462   hCos2PhiBkgLoScal->SetLineStyle(7);
 463   hCos2PhiBkgHiScal->SetLineStyle(2);
 464   hCos2PhiBkgLoScal->SetLineColor(kRed+1);
 465   hCos2PhiBkgHiScal->SetLineColor(kBlue+1);
 466   printf("After Scaling %d: BkgLo = %f+-%f BkgHi=%f+-%f \n",theBinR,
 467          hCos2PhiBkgLoScal->GetBinContent(theBinR),hCos2PhiBkgLoScal->GetBinError(theBinR),
 468          hCos2PhiBkgHiScal->GetBinContent(theBinR),hCos2PhiBkgHiScal->GetBinError(theBinR));
 469
 470
 471   TH1F* hCos2PhiBkgAver=0x0;
 472   if(optBkg==0){
 473     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 474     hCos2PhiBkgAver->Add(hCos2PhiBkgHiScal);
 475     hCos2PhiBkgAver->Scale(0.5);
 476   }else if(optBkg==-1){
 477     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 478   }else{
 479     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgHiScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 480   }
 481   printf("After Average %d: BkgLo = %f+-%f BkgHi=%f+-%f BkgAve=%f+-%f \n",theBinR,
 482          hCos2PhiBkgLoScal->GetBinContent(theBinR),hCos2PhiBkgLoScal->GetBinError(theBinR),
 483          hCos2PhiBkgHiScal->GetBinContent(theBinR),hCos2PhiBkgHiScal->GetBinError(theBinR),
 484          hCos2PhiBkgAver->GetBinContent(theBinR),hCos2PhiBkgAver->GetBinError(theBinR)
 485          );
 486
 487   hCos2PhiBkgAver->SetLineWidth(2);
 488   hCos2PhiBkgAver->SetLineColor(kGreen+1);
 489   hCos2PhiBkgAver->SetFillColor(kGreen+1);
 490   TH1F* hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone(Form("hCos2PhiSigBin%d",iFinalPtBin));
 491   hCos2PhiSig->Add(hCos2PhiBkgAver,-1.);
 492   printf("After Subtraction %d: All = %f+-%f BkgAve=%f+-%f Sig=%f+-%f \n",theBinR,
 493          hCos2PhiSigReg->GetBinContent(theBinR),hCos2PhiSigReg->GetBinError(theBinR),
 494          hCos2PhiBkgAver->GetBinContent(theBinR),hCos2PhiBkgAver->GetBinError(theBinR),
 495          hCos2PhiSig->GetBinContent(theBinR),hCos2PhiSig->GetBinError(theBinR)
 496          );
 497
 498   TLegend* leg0=new TLegend(0.3,0.6,0.75,0.89);
 499   TPaveText* t0= new TPaveText(0.15,0.70,0.45,0.89,"NDC");
 500   if(c1){
 501     c1->cd(3);
 502     hCos2PhiSigReg->Draw();
 503     hCos2PhiBkgAver->Draw("same");
 504     hCos2PhiBkgLoScal->Draw("same");
 505     hCos2PhiBkgHiScal->Draw("same");
 506     leg0->SetFillColor(0);
 507     TLegendEntry* ent=leg0->AddEntry(hCos2PhiSigReg,"Signal region","L");
 508     ent->SetTextColor(hCos2PhiSigReg->GetLineColor());
 509     ent=leg0->AddEntry(hCos2PhiBkgLoScal,"Left side band","L");
 510     ent->SetTextColor(hCos2PhiBkgLoScal->GetLineColor());
 511     ent=leg0->AddEntry(hCos2PhiBkgHiScal,"Right side band","L");
 512     ent->SetTextColor(hCos2PhiBkgHiScal->GetLineColor());
 513     ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","F");
 514     ent->SetTextColor(hCos2PhiBkgAver->GetLineColor());
 515     leg0->Draw();
 516     c1->cd(4);
 517     hCos2PhiSig->Draw("EP");
 518     t0->SetFillColor(0);
 519     t0->AddText(Form("v2=%.3f+-%.3f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError()));
 520     t0->Draw();
 521
 522     c1->Update();
 523   }
 524   if(!c1){
 525     delete leg0;
 526     delete t0;
 527     delete hCos2PhiBkgLo;
 528     delete hCos2PhiBkgHi;
 529     delete hCos2PhiSigReg;
 530     delete hCos2PhiBkgLoScal;
 531     delete hCos2PhiBkgHiScal;
 532     delete hCos2PhiBkgAver;
 533   }
 534   printf("Signal from mass fitter = %f  Signal from subracted histo= %f\n",
 535          signalFromFit,hCos2PhiSig->Integral());
 536   return hCos2PhiSig;
 537 }
 538
 539
 540 //______________________________________________________________________________
 541 TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors, Bool_t useConst){
 542
 543   Int_t npars=fSB->GetNpar();
 544   Double_t sigmaSB=fSB->GetParameter(npars-1);
 545   Double_t massSB=fSB->GetParameter(npars-2);
 546   Double_t integr=fSB->GetParameter(npars-3);
 547   Double_t sOverAll=(fSB->Eval(massSB)-fB2->Eval(massSB))/fSB->Eval(massSB);
 548   printf("mass=%f  S+B=%f   bkg=%f S/(S+B)=%f\n",massSB,fSB->Eval(massSB),fB2->Eval(massSB),sOverAll);
 549   printf("Number of parameters: %d. Signal params: %f %f %f\n",npars,massSB,sigmaSB,integr);
 550   if(optErrors==1){
 551     sOverAll+=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
 552   }
 553   else if(optErrors==-1){
 554     sOverAll-=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
 555   }
 556
 557   Int_t nbinsmass=hMassDphi->GetNbinsY();
 558   Double_t minmass=hMassDphi->GetYaxis()->GetXmin();
 559   Double_t maxmass=hMassDphi->GetYaxis()->GetXmax();
 560
 561   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);
 562   fSig->SetParameters(integr,massSB,sigmaSB);
 563
 564   TH1F* hAverCos2Phi=new TH1F(Form("hAverCos2PhiBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 565   TH1F* hFractionSig=new TH1F(Form("hFractionSigBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 566   TH1F* hFractionBkg=new TH1F(Form("hFractionBkgBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 567   for(Int_t iBin=1; iBin<=hMassDphi->GetNbinsY(); iBin++){
 568     TH1F* htemp=(TH1F*)hMassDphi->ProjectionX("htemp",iBin,iBin);
 569     hAverCos2Phi->SetBinContent(iBin,htemp->GetMean());
 570     hAverCos2Phi->SetBinError(iBin,htemp->GetMeanError());
 571     Double_t sig=fSig->Eval(hFractionSig->GetBinCenter(iBin));
 572     Double_t bkg=fB2->Eval(hFractionSig->GetBinCenter(iBin));
 573     if(bkg<1 && sig<1){
 574       hFractionSig->SetBinContent(iBin,0.);
 575       hFractionSig->SetBinError(iBin,0.);
 576       hFractionBkg->SetBinContent(iBin,1.);
 577       hFractionBkg->SetBinError(iBin,0.);
 578     }else{
 579       Double_t fracs=sig/(sig+bkg);
 580       Double_t fracb=bkg/(sig+bkg);
 581       Double_t efracs=0.;//TMath::Sqrt(fracs*(1.-fracs)/(sig+bkg));
 582       Double_t efracb=0.;//TMath::Sqrt(fracb*(1.-fracb)/(sig+bkg));
 583
 584       hFractionSig->SetBinContent(iBin,fracs);
 585       hFractionSig->SetBinError(iBin,efracs);
 586       hFractionBkg->SetBinContent(iBin,fracb);
 587       hFractionBkg->SetBinError(iBin,efracb);
 588     }
 589     delete htemp;
 590   }
 591
 592   TF1* fv2=0x0;
 593   if(useConst){
 594     fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMass,minmass,maxmass,5);
 595   }else{
 596     fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMassLin,minmass,maxmass,6);
 597     fv2->SetParameter(5,0.);
 598   }
 599   fv2->SetParameter(0,sOverAll);
 600   fv2->SetParameter(1,massSB);
 601   fv2->SetParameter(2,sigmaSB);
 602   fv2->SetParameter(3,0.2);
 603   fv2->SetParameter(4,0.2);
 604   fv2->FixParameter(0,sOverAll);
 605   fv2->FixParameter(1,massSB);
 606   fv2->FixParameter(2,sigmaSB);
 607
 608   printf("LOOP %d\n",rebin);
 609   if((hAverCos2Phi->GetNbinsX()%rebin)==0){
 610     hAverCos2Phi->Rebin(rebin);
 611     hAverCos2Phi->Scale(1./(Double_t)rebin);
 612   }
 613
 614   if(c2){
 615     c2->Divide(2,2);
 616     c2->cd(1);
 617     hMassDphi->Draw("colz");
 618     c2->cd(2);
 619     hMass->Rebin(2);
 620     hMass->SetMinimum(0.);
 621     hMass->SetMarkerStyle(20);
 622     hMass->Draw("E");
 623     fSB->Draw("same");
 624     fSig->Draw("same");
 625     fB2->Draw("same");
 626     c2->cd(3);
 627     hFractionSig->SetMaximum(1.2);
 628     hFractionSig->Draw();
 629     hFractionSig->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 630     hFractionSig->GetYaxis()->SetTitle("Fraction");
 631     hFractionBkg->SetLineColor(2);
 632     hFractionBkg->Draw("same");
 633     TLegend* leg1=new TLegend(0.15,0.15,0.35,0.35);
 634     leg1->SetFillColor(0);
 635     TLegendEntry* ent=leg1->AddEntry(hFractionSig,"S/(S+B)","L");
 636     ent->SetTextColor(hFractionSig->GetLineColor());
 637     ent=leg1->AddEntry(hFractionBkg,"B/(S+B)","L");
 638     ent->SetTextColor(hFractionBkg->GetLineColor());
 639     leg1->Draw();
 640     c2->cd(4);
 641     hAverCos2Phi->Fit(fv2);
 642     fv2->DrawCopy("same");
 643     hAverCos2Phi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 644     hAverCos2Phi->GetYaxis()->SetTitle("v_2^{obs}");
 645     TPaveText* t1= new TPaveText(0.55,0.70,0.89,0.89,"NDC");
 646     t1->SetFillColor(0);
 647     t1->AddText(Form("v2sig=%.3f+-%.3f\n",fv2->GetParameter(3),fv2->GetParError(3)));
 648     if(useConst){
 649       t1->AddText(Form("v2bkg=%.3f+-%.3f\n",fv2->GetParameter(4),fv2->GetParError(4)));
 650     }else{
 651       t1->AddText(Form("v2bkg=(%.3f+-%.3f) + (%.3g+-%.3g)*mass\n",fv2->GetParameter(4),fv2->GetParError(4),fv2->GetParameter(5),fv2->GetParError(5)));
 652     }
 653     t1->Draw();
 654     c2->Update();
 655   }else{
 656     hAverCos2Phi->Fit(fv2);
 657   }
 658   return fv2;
 659 }
 660
 661
 662 //______________________________________________________________________________
 663 Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh){
 664   // Event plane resolution syst err (from wide centrality bin
 665   TH1F* hResolSubAB=0x0;
 666   TH1F* hResolSubBC=0x0;
 667   TH1F* hResolSubAC=0x0;
 668   Double_t xmin=1.;
 669   Double_t xmax=-1.;
 670   TGraphAsymmErrors* grSingle=new TGraphAsymmErrors(0);
 671   TGraphAsymmErrors* grInteg=new TGraphAsymmErrors(0);
 672   Int_t iPt=0;
 673   Int_t minCentTimesTen=minCent*10;
 674   Int_t maxCentTimesTen=maxCent*10;
 675   Double_t binHalfWid=25./20.;
 676   for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){
 677     TString hisnameAB=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+25);
 678     TString hisnameBC=Form("hEvPlaneReso2centr%d_%d",iHisC,iHisC+25);
 679     TString hisnameAC=Form("hEvPlaneReso3centr%d_%d",iHisC,iHisC+25);
 680     TH1F* hResolSubABsing=(TH1F*)lst->FindObject(hisnameAB.Data());
 681     cout<<hResolSubABsing->GetName()<<endl;
 682     TH1F* hResolSubBCsing=0x0;
 683     TH1F* hResolSubACsing=0x0;
 684     if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 685        evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
 686       hResolSubBCsing=(TH1F*)lst->FindObject(hisnameBC.Data());
 687       cout<<hResolSubBCsing->GetName()<<endl;
 688       hResolSubACsing=(TH1F*)lst->FindObject(hisnameAC.Data());
 689       cout<<hResolSubACsing->GetName()<<endl;
 690     }
 691     Double_t error;
 692     Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubABsing,hResolSubBCsing,hResolSubACsing);
 693     Double_t resolFullmin=resolFull-error;
 694     Double_t resolFullmax=resolFull+error;
 695
 696     Double_t binCentr=(Double_t)iHisC/10.+binHalfWid;
 697     grSingle->SetPoint(iPt,binCentr,resolFull);
 698     grSingle->SetPointEXlow(iPt,binHalfWid);
 699     grSingle->SetPointEXhigh(iPt,binHalfWid);
 700     grSingle->SetPointEYlow(iPt,resolFullmax-resolFull);
 701     grSingle->SetPointEYhigh(iPt,resolFull-resolFullmin);
 702     ++iPt;
 703     if(resolFullmin<xmin) xmin=resolFullmin;
 704     if(resolFullmax>xmax) xmax=resolFullmax;
 705     if(iHisC==minCentTimesTen){
 706       hResolSubAB=(TH1F*)hResolSubABsing->Clone("hResolSubAB");
 707       if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 708          evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
 709         hResolSubBC=(TH1F*)hResolSubBCsing->Clone("hResolSubAB");
 710         hResolSubAC=(TH1F*)hResolSubACsing->Clone("hResolSubAB");
 711       }
 712     }else{
 713       hResolSubAB->Add(hResolSubABsing);
 714       if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 715          evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
 716         hResolSubBC->Add(hResolSubBCsing);
 717         hResolSubAC->Add(hResolSubACsing);
 718       }
 719     }
 720     printf("Adding histogram %s\n",hResolSubAB->GetName());
 721   }
 722   rcflow=xmin;
 723   rcfhigh=xmax;
 724
 725
 726   Double_t error;
 727   Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubAB,hResolSubBC,hResolSubAC);
 728
 729   grInteg->SetPoint(0,40.,resolFull);
 730   grInteg->SetPointEXlow(0,10);
 731   grInteg->SetPointEXhigh(0,10.);
 732   grInteg->SetPointEYlow(0,error);
 733   grInteg->SetPointEYhigh(0,error);
 734
 735   TCanvas* cEP=new TCanvas("cEP","EvPlaneResol");
 736   cEP->Divide(1,2);
 737   cEP->cd(1);
 738   hResolSubAB->Draw();
 739   if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
 740      evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
 741     hResolSubBC->SetLineColor(2);
 742     hResolSubBC->Draw("same");
 743     hResolSubAC->SetLineColor(4);
 744     hResolSubAC->Draw("same");
 745   }
 746   TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
 747   tres->SetNDC();
 748   tres->Draw();
 749   cEP->cd(2);
 750   grSingle->SetMarkerStyle(20);
 751   grInteg->SetMarkerColor(kRed+1);
 752   grInteg->SetLineColor(kRed+1);
 753   grInteg->SetMarkerStyle(29);
 754   grSingle->Draw("AP");
 755   grSingle->GetXaxis()->SetTitle("Centrality Percentile");
 756   grSingle->GetYaxis()->SetTitle("Resolution Correction Factor");
 757   grInteg->Draw("Psame");
 758
 759   return resolFull;
 760
 761 }
 762
 763 //______________________________________________________________________________
 764 Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
 765   Double_t resolFull;
 766   if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub ||
 767      evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){
 768     resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
 769     error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
 770   }else if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
 771     if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta){
 772       resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
 773       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
 774     }else if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
 775       resolFull=AliVertexingHFUtils::GetSubEvResol(hsubev1);
 776       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hsubev1));
 777     }
 778   }else{
 779     Double_t resolSub[3];
 780     Double_t errors[3];
 781     TH1F* hevplresos[3];
 782     hevplresos[0]=hsubev1;
 783     hevplresos[1]=hsubev2;
 784     hevplresos[2]=hsubev3;
 785     for(Int_t ires=0;ires<3;ires++){
 786       resolSub[ires]=hevplresos[ires]->GetMean();
 787       errors[ires]=hevplresos[ires]->GetMeanError();
 788     }
 789     Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
 790     if(evPlane==AliEventPlaneResolutionHandler::kVZEROC ||
 791        evPlane==AliEventPlaneResolutionHandler::kVZERO){
 792       resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
 793       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
 794     }
 795     else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){
 796       resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
 797       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
 798     }
 799     else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta ||
 800             evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
 801       resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
 802       error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
 803     }
 804   }
 805   return resolFull;
 806 }
 807
 808 //______________________________________________________________________________
 809 void LoadMassHistos(TList* lst, TH2F** hMassDphi){
 810
 811   Int_t minCentTimesTen=minCent*10;
 812   Int_t maxCentTimesTen=maxCent*10;
 813   for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){
 814     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 815       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){
 816         TString hisname=Form("hMc2deltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
 817         TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
 818         if(hMassDphi[iFinalPtBin]==0x0){
 819           hMassDphi[iFinalPtBin]=(TH2F*)htmp->Clone(Form("hMassCos2DphiBin%d",iFinalPtBin));
 820         }else{
 821           hMassDphi[iFinalPtBin]->Add(htmp);
 822         }
 823         printf("Adding histogram %s to PtBin %d\n",hisname.Data(),iFinalPtBin);
 824       }
 825     }
 826   }
 827 }
 828
 829 //______________________________________________________________________________
 830 Bool_t DefinePtBins(TDirectoryFile* df){
 831
 832   AliRDHFCuts *d0cuts;
 833   if(partname.Contains("Dzero")) d0cuts=(AliRDHFCutsD0toKpi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
 834   if(partname.Contains("Dplus")) d0cuts=(AliRDHFCutsDplustoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
 835   if(partname.Contains("Dstar")) d0cuts=((AliRDHFCutsDStartoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName()));
 836
 837   Int_t nPtBinsCuts=d0cuts->GetNPtBins();
 838   printf("Number of pt bins for cut object = %d\n",nPtBinsCuts);
 839   Float_t *ptlimsCuts=d0cuts->GetPtBinLimits();
 840   for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
 841   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
 842     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 843       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[iFinalPtBin])<0.0001){
 844         minPtBin[iFinalPtBin]=iPtCuts;
 845         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
 846       }
 847     }
 848     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[nFinalPtBins])<0.0001) maxPtBin[nFinalPtBins-1]=iPtCuts-1;
 849   }
 850   if(TMath::Abs(ptlims[nFinalPtBins]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nFinalPtBins-1]=nPtBinsCuts-1;
 851   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 852     printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
 853     if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
 854   }
 855
 856   return kTRUE;
 857 }
 858
 859
 860 //______________________________________________________________________________
 861 void SystForSideBands(){
 862   // Compute systematic ucnertainty for the side band subtraction method
 863   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
 864
 865   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
 866   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
 867
 868   TFile *f = TFile::Open(filename.Data());
 869   if(!f){
 870     printf("file %s not found, please check file name\n",filename.Data());
 871     return;
 872   }
 873   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 874   if(!dir){
 875     printf("Directory %s not found, please check dir name\n",dirname.Data());
 876     return;
 877   }
 878   Bool_t binOK=DefinePtBins(dir);
 879   if(!binOK){
 880     printf("ERROR: mismatch in pt binning\n");
 881     return;
 882   }
 883   TList *lst =(TList*)dir->Get(listname.Data());
 884   if(!lst){
 885     printf("list %s not found in file, please check list name\n",listname.Data());
 886     return;
 887   }
 888
 889   AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
 890   epres->SetEventPlane(evPlane);
 891   epres->SetResolutionOption(evPlaneRes);
 892   epres->SetUseNcollWeights(useNcollWeight);
 893   Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
 894   delete epres;
 895   printf("Event plane resolution %f\n",resolFull);
 896
 897   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
 898   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 899     hMassDphi[iFinalPtBin]=0x0;
 900   }
 901   LoadMassHistos(lst, hMassDphi);
 902
 903   Int_t nSteps=21;
 904
 905   TGraphErrors** gSystSigRange=new TGraphErrors*[nFinalPtBins];
 906   TGraphErrors** gSystBkgRange=new TGraphErrors*[nFinalPtBins];
 907   TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
 908   TGraphErrors** gSystWhichSide=new TGraphErrors*[nFinalPtBins];
 909
 910   Double_t min1[nFinalPtBins],max1[nFinalPtBins];
 911   Double_t min123[nFinalPtBins],max123[nFinalPtBins];
 912   Double_t min2[nFinalPtBins],max2[nFinalPtBins];
 913   Double_t min3[nFinalPtBins],max3[nFinalPtBins];
 914   Double_t min4[nFinalPtBins],max4[nFinalPtBins];
 915
 916   Double_t sigmaRangeForSigOrig=sigmaRangeForSig;
 917   Double_t sigmaRangeForBkgOrig=sigmaRangeForBkg;
 918
 919
 920   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 921     printf("**************** BIN %d ******************\n",iFinalPtBin);
 922
 923     Int_t rebinHistoSideBandsOrig=rebinHistoSideBands[iFinalPtBin];
 924
 925     gSystSigRange[iFinalPtBin]=new TGraphErrors(0);
 926     gSystSigRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Signal Region Ptbin %d",iFinalPtBin));
 927
 928     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 929
 930     Bool_t out=FitMassSpectrum(hMass,0x0);
 931     if(!out) continue;
 932
 933     min1[iFinalPtBin]=99999.;
 934     max1[iFinalPtBin]=-99999.;
 935     min123[iFinalPtBin]=99999.;
 936     max123[iFinalPtBin]=-99999.;
 937     for(Int_t iStep=0; iStep<nSteps; iStep++){
 938       Int_t index=iFinalPtBin*nSteps+iStep;
 939       sigmaRangeForSig=1.5+0.1*iStep;
 940       sigmaRangeForBkg=sigmaRangeForBkgOrig;
 941       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 942       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 943       Double_t v2obsM1=hCos2PhiSig->GetMean();
 944       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 945       delete hCos2PhiSig;
 946       Double_t v2M1=v2obsM1/resolFull;
 947       Double_t errv2M1=errv2obsM1/resolFull;
 948       gSystSigRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForSig,v2M1);
 949       gSystSigRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
 950       if(v2M1>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M1;
 951       if(v2M1<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M1;
 952       if(sigmaRangeForSig>=2. && sigmaRangeForSig<=3){
 953         if(v2M1>max123[iFinalPtBin]) max123[iFinalPtBin]=v2M1;
 954         if(v2M1<min123[iFinalPtBin]) min123[iFinalPtBin]=v2M1;
 955       }
 956     }
 957
 958     min2[iFinalPtBin]=99999.;
 959     max2[iFinalPtBin]=-99999.;
 960     gSystBkgRange[iFinalPtBin]=new TGraphErrors(0);
 961     gSystBkgRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Bkg Region Ptbin %d",iFinalPtBin));
 962     for(Int_t iStep=0; iStep<nSteps; iStep++){
 963       Int_t index=nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
 964       sigmaRangeForSig=sigmaRangeForSigOrig;
 965       sigmaRangeForBkg=4.+0.1*iStep;
 966       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
 967       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 968       Double_t v2obsM1=hCos2PhiSig->GetMean();
 969       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 970       delete hCos2PhiSig;
 971       Double_t v2M1=v2obsM1/resolFull;
 972       Double_t errv2M1=errv2obsM1/resolFull;
 973       gSystBkgRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForBkg,v2M1);
 974       gSystBkgRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
 975       if(v2M1>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M1;
 976       if(v2M1<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M1;
 977     }
 978
 979     min3[iFinalPtBin]=99999.;
 980     max3[iFinalPtBin]=-99999.;
 981     gSystRebin[iFinalPtBin]=new TGraphErrors(0);
 982     gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin Ptbin %d",iFinalPtBin));
 983     Int_t nPts=0;
 984     for(Int_t iStep=0; iStep<nSteps; iStep++){
 985       Int_t index=2*nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
 986       sigmaRangeForSig=sigmaRangeForSigOrig;
 987       sigmaRangeForBkg=sigmaRangeForBkgOrig;
 988       rebinHistoSideBands[iFinalPtBin]=iStep+1;
 989       if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistoSideBands[iFinalPtBin]!=0) continue;
 990       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
 991       Double_t v2obsM1=hCos2PhiSig->GetMean();
 992       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 993       delete hCos2PhiSig;
 994       Double_t v2M1=v2obsM1/resolFull;
 995       Double_t errv2M1=errv2obsM1/resolFull;
 996       gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistoSideBands[iFinalPtBin],v2M1);
 997       gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
 998       if(v2M1>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M1;
 999       if(v2M1<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M1;
1000       ++nPts;
1001     }
1002
1003     min4[iFinalPtBin]=99999.;
1004     max4[iFinalPtBin]=-99999.;
1005     gSystWhichSide[iFinalPtBin]=new TGraphErrors(0);
1006     gSystWhichSide[iFinalPtBin]->SetTitle(Form("v2 vs. WhichSide Ptbin %d",iFinalPtBin));
1007     nPts=0;
1008     for(Int_t iStep=-1; iStep<=1; iStep++){
1009       Int_t index=3*nSteps*nFinalPtBins+2+iStep;
1010       sigmaRangeForSig=sigmaRangeForSigOrig;
1011       sigmaRangeForBkg=sigmaRangeForBkgOrig;
1012       rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
1013       TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0,iStep);
1014       Double_t v2obsM1=hCos2PhiSig->GetMean();
1015       Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
1016       delete hCos2PhiSig;
1017       Double_t v2M1=v2obsM1/resolFull;
1018       Double_t errv2M1=errv2obsM1/resolFull;
1019       gSystWhichSide[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M1);
1020       gSystWhichSide[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
1021       if(v2M1>max4[iFinalPtBin]) max4[iFinalPtBin]=v2M1;
1022       if(v2M1<min4[iFinalPtBin]) min4[iFinalPtBin]=v2M1;
1023       ++nPts;
1024     }
1025     rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
1026   }
1027
1028   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1029     printf("------ Pt Bin %d ------\n",iFinalPtBin);
1030     printf("Range of values for sig variation = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
1031     printf("           (limited to 2-3 sigma) = %f %f\n",min123[iFinalPtBin],max123[iFinalPtBin]);
1032     printf("Range of values for bkg variation = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
1033     printf("Range of values for rebin = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
1034     printf("Range of values for whichside = %f %f\n",min4[iFinalPtBin],max4[iFinalPtBin]);
1035     Float_t minenv=TMath::Min(min123[iFinalPtBin],TMath::Min(min2[iFinalPtBin],min4[iFinalPtBin]));
1036     Float_t maxenv=TMath::Max(max123[iFinalPtBin],TMath::Max(max2[iFinalPtBin],max4[iFinalPtBin]));
1037     printf(" ---> Envelope                = %f %f\n",minenv,maxenv);
1038   }
1039
1040   TCanvas* cs1=new TCanvas("cs1");
1041   cs1->Divide(2,2);
1042   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1043     cs1->cd(iFinalPtBin+1);
1044     gSystSigRange[iFinalPtBin]->SetMarkerStyle(20);
1045     gSystSigRange[iFinalPtBin]->Draw("AP");
1046     gSystSigRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaSignal");
1047     gSystSigRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1048   }
1049
1050   TCanvas* cs2=new TCanvas("cs2");
1051   cs2->Divide(2,2);
1052   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1053     cs2->cd(iFinalPtBin+1);
1054     gSystBkgRange[iFinalPtBin]->SetMarkerStyle(20);
1055     gSystBkgRange[iFinalPtBin]->Draw("AP");
1056     gSystBkgRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaBackground");
1057     gSystBkgRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1058   }
1059
1060   TCanvas* cs3=new TCanvas("cs3");
1061   cs3->Divide(2,2);
1062   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1063     cs3->cd(iFinalPtBin+1);
1064     gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
1065     gSystRebin[iFinalPtBin]->Draw("AP");
1066     gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
1067     gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1068   }
1069
1070   TCanvas* cs4=new TCanvas("cs4");
1071   cs4->Divide(2,2);
1072   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1073     cs4->cd(iFinalPtBin+1);
1074     gSystWhichSide[iFinalPtBin]->SetMarkerStyle(20);
1075     gSystWhichSide[iFinalPtBin]->Draw("AP");
1076     gSystWhichSide[iFinalPtBin]->GetXaxis()->SetTitle("Side band used");
1077     gSystWhichSide[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1078   }
1079 }
1080
1081 //______________________________________________________________________________
1082 void SystForFitv2Mass(){
1083   // Compute systematic ucnertainty for the v2(M) method
1084   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
1085
1086   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
1087   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
1088
1089   TFile *f = TFile::Open(filename.Data());
1090   if(!f){
1091     printf("file %s not found, please check file name\n",filename.Data());
1092     return;
1093   }
1094   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
1095   if(!dir){
1096     printf("Directory %s not found, please check dir name\n",dirname.Data());
1097     return;
1098   }
1099   Bool_t binOK=DefinePtBins(dir);
1100   if(!binOK){
1101     printf("ERROR: mismatch in pt binning\n");
1102     return;
1103   }
1104   TList *lst =(TList*)dir->Get(listname.Data());
1105   if(!lst){
1106     printf("list %s not found in file, please check list name\n",listname.Data());
1107     return;
1108   }
1109
1110   AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
1111   epres->SetEventPlane(evPlane);
1112   epres->SetResolutionOption(evPlaneRes);
1113   epres->SetUseNcollWeights(useNcollWeight);
1114   Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
1115   delete epres;
1116   printf("Event plane resolution %f\n",resolFull);
1117
1118   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
1119   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1120     hMassDphi[iFinalPtBin]=0x0;
1121   }
1122   LoadMassHistos(lst, hMassDphi);
1123
1124   Int_t nSteps=11;
1125
1126   TGraphErrors** gSystParamErr=new TGraphErrors*[nFinalPtBins];
1127   TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
1128   TGraphErrors** gSystLinConst=new TGraphErrors*[nFinalPtBins];
1129
1130   Double_t min1[nFinalPtBins],max1[nFinalPtBins];
1131   Double_t min2[nFinalPtBins],max2[nFinalPtBins];
1132   Double_t min3[nFinalPtBins],max3[nFinalPtBins];
1133
1134
1135   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1136     printf("**************** BIN %d ******************\n",iFinalPtBin);
1137
1138     Int_t rebinHistov2MassOrig=rebinHistov2Mass[iFinalPtBin];
1139     TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
1140
1141     Bool_t out=FitMassSpectrum(hMass,0x0);
1142     if(!out) continue;
1143
1144     min1[iFinalPtBin]=99999.;
1145     max1[iFinalPtBin]=-99999.;
1146     gSystRebin[iFinalPtBin]=new TGraphErrors(0);
1147     gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin PtBin %d",iFinalPtBin));
1148     Int_t nPts=0;
1149     for(Int_t iStep=0; iStep<nSteps; iStep++){
1150       Int_t index=iStep;
1151       rebinHistov2Mass[iFinalPtBin]=iStep+1;
1152       if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistov2Mass[iFinalPtBin]!=0) continue;
1153       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,useConstantvsBkgVsMass);
1154       Double_t v2obsM2=fv2->GetParameter(3);
1155       Double_t errv2obsM2=fv2->GetParError(3);
1156       delete fv2;
1157       Double_t v2M2=v2obsM2/resolFull;
1158       Double_t errv2M2=errv2obsM2/resolFull;
1159       gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistov2Mass[iFinalPtBin],v2M2);
1160       gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
1161       if(v2M2>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M2;
1162       if(v2M2<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M2;
1163       ++nPts;
1164     }
1165     rebinHistov2Mass[iFinalPtBin]=rebinHistov2MassOrig;
1166
1167     min2[iFinalPtBin]=99999.;
1168     max2[iFinalPtBin]=-99999.;
1169     gSystParamErr[iFinalPtBin]=new TGraphErrors(0);
1170     gSystParamErr[iFinalPtBin]->SetTitle(Form("v2 vs. ParamErr PtBin %d",iFinalPtBin));
1171     nPts=0;
1172     for(Int_t iStep=-1; iStep<=1; iStep++){
1173       Int_t index=nSteps*2+iStep;
1174       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,iStep,useConstantvsBkgVsMass);
1175       Double_t v2obsM2=fv2->GetParameter(3);
1176       Double_t errv2obsM2=fv2->GetParError(3);
1177       delete fv2;
1178       Double_t v2M2=v2obsM2/resolFull;
1179       Double_t errv2M2=errv2obsM2/resolFull;
1180       gSystParamErr[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M2);
1181       gSystParamErr[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
1182       if(v2M2>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M2;
1183       if(v2M2<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M2;
1184       ++nPts;
1185     }
1186
1187     min3[iFinalPtBin]=99999.;
1188     max3[iFinalPtBin]=-99999.;
1189     gSystLinConst[iFinalPtBin]=new TGraphErrors(0);
1190     gSystLinConst[iFinalPtBin]->SetTitle(Form("v2 LinVsConst Bin%d",iFinalPtBin));
1191     nPts=0;
1192     for(Int_t iStep=0; iStep<=1; iStep++){
1193       Int_t index=nSteps*3+iStep;
1194       TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,iStep);
1195       Double_t v2obsM2=fv2->GetParameter(3);
1196       Double_t errv2obsM2=fv2->GetParError(3);
1197       delete fv2;
1198       Double_t v2M2=v2obsM2/resolFull;
1199       Double_t errv2M2=errv2obsM2/resolFull;
1200       gSystLinConst[iFinalPtBin]->SetPoint(nPts,1.-(Double_t)iStep,v2M2);
1201       gSystLinConst[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
1202       if(v2M2>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M2;
1203       if(v2M2<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M2;
1204       ++nPts;
1205     }
1206
1207
1208   }
1209
1210   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1211     printf("------ Pt Bin %d ------\n",iFinalPtBin);
1212     printf("Range of values for rebin = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
1213     printf("Range of values for par err = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
1214     printf("Range of values for lin const = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
1215     Float_t minenv=TMath::Min(min2[iFinalPtBin],min3[iFinalPtBin]);
1216     Float_t maxenv=TMath::Max(max2[iFinalPtBin],max3[iFinalPtBin]);
1217     printf(" ---> Envelope                = %f %f\n",minenv,maxenv);
1218   }
1219
1220   TCanvas* cs1=new TCanvas("cs1");
1221   cs1->Divide(2,2);
1222   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1223     cs1->cd(iFinalPtBin+1);
1224     gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
1225     gSystRebin[iFinalPtBin]->Draw("AP");
1226     gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
1227     gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1228   }
1229
1230   TCanvas* cs2=new TCanvas("cs2");
1231   cs2->Divide(2,2);
1232   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1233     cs2->cd(iFinalPtBin+1);
1234     gSystParamErr[iFinalPtBin]->SetMarkerStyle(20);
1235     gSystParamErr[iFinalPtBin]->Draw("AP");
1236     gSystParamErr[iFinalPtBin]->GetXaxis()->SetTitle("Error on Signal yield");
1237     gSystParamErr[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1238   }
1239
1240   TCanvas* cs3=new TCanvas("cs3");
1241   cs3->Divide(2,2);
1242   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
1243     cs3->cd(iFinalPtBin+1);
1244     gSystLinConst[iFinalPtBin]->SetMarkerStyle(20);
1245     gSystLinConst[iFinalPtBin]->Draw("AP");
1246     gSystLinConst[iFinalPtBin]->GetXaxis()->SetLimits(-0.1,1.1);
1247     gSystLinConst[iFinalPtBin]->GetXaxis()->SetTitle("Const/Linear v2 of background");
1248     gSystLinConst[iFinalPtBin]->GetYaxis()->SetTitle("v2");
1249   }
1250 }