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
  20 #include "AliAODRecoDecayHF.h"
  21 #include "AliRDHFCuts.h"
  22 #include "AliRDHFCutsDplustoKpipi.h"
  23 #include "AliRDHFCutsD0toKpi.h"
  24 #include "AliHFMassFitter.h"
  25 #include "AliEventPlaneResolution.h"
  26 #endif
  27
  28 Double_t v2vsMass(Double_t *x, Double_t *par){
  29   // Fit function for signal+background
  30   // par[0] = S/B at the mass peak
  31   // par[1] = mass
  32   // par[2] = sigma
  33   // par[3] = v2sig
  34   // par[4] = v2back
  35
  36   Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
  37   Double_t fracbkg=1-fracsig;
  38   return par[3]*fracsig+par[4]*fracbkg;
  39 }
  40
  41 void Extractv2from2Dhistos(){
  42
  43   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
  44
  45   TFile* fil=new TFile("AnalysisResults_ptbins.root");
  46   TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
  47   TList* lst=(TList*)df->Get("coutputv2D0Std");
  48
  49   AliRDHFCutsD0toKpi *d0cuts=(AliRDHFCutsD0toKpi*)df->Get("D0toKpiCuts");
  50   Int_t nPtBinsCuts=d0cuts->GetNPtBins();
  51   printf("Number of pt bins for cut object = %d\n",nPtBinsCuts);
  52   Float_t *ptlimsCuts=d0cuts->GetPtBinLimits();
  53   for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f \n",iPt,ptlimsCuts[iPt]);
  54
  55   Int_t minCent=30;
  56   Int_t maxCent=50;
  57
  58
  59   const Int_t nFinalPtBins=3;
  60   Double_t ptlims[nFinalPtBins+1]={2.,4.,5.,12.};
  61
  62   Int_t minPtBin[nFinalPtBins]={-1,-1,-1};
  63   Int_t maxPtBin[nFinalPtBins]={-1,-1,-1};
  64   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
  65     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
  66       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[iFinalPtBin])<0.0001){
  67         minPtBin[iFinalPtBin]=iPtCuts;
  68         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts;
  69       }
  70     }
  71     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[nFinalPtBins])<0.0001) maxPtBin[nFinalPtBins-1]=iPtCuts;
  72  }
  73   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++) printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
  74
  75   TH1F* hResolSubAB=0x0;
  76   TH2F** hMassDphi=new TH2F*[nFinalPtBins];
  77   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
  78     hMassDphi[iFinalPtBin]=0x0;
  79   }
  80
  81
  82   for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){
  83     TString hisnameEP=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+5);
  84     TH2F* htmpEP=(TH2F*)lst->FindObject(hisnameEP.Data());
  85
  86     if(iHisC==minCent){
  87       hResolSubAB=(TH1F*)htmpEP->Clone("hResolSubAB");
  88     }else{
  89       hResolSubAB->Add(htmpEP);
  90     }
  91     printf("Adding histogram %s\n",hisnameEP.Data());
  92
  93     for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
  94       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){
  95         TString hisname=Form("hMc2phi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+5);
  96         TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
  97         if(hMassDphi[iFinalPtBin]==0x0){
  98           hMassDphi[iFinalPtBin]=(TH2F*)htmp->Clone(Form("hMassCos2DphiBin%d",iFinalPtBin));
  99         }else{
 100           hMassDphi[iFinalPtBin]->Add(htmp);
 101         }
 102         printf("Adding histogram %s to PtBin %d\n",hisname.Data(),iFinalPtBin);
 103
 104       }
 105     }
 106   }
 107
 108   AliEventPlaneResolution *resol=new AliEventPlaneResolution(1);
 109   resol->SetSubEventHisto(hResolSubAB);
 110
 111   Double_t resolSub=AliEventPlaneResolution::GetSubEvResol(hResolSubAB);
 112   printf("Resolution on sub event  = %.4f\n",resolSub);
 113   Double_t chisub=AliEventPlaneResolution::FindChi(resolSub,1);
 114   printf("Chi Subevent             = %.4f\n",chisub);
 115   Double_t chifull=chisub*TMath::Sqrt(2);
 116   printf("Chi Full Event           = %.4f\n",chifull);
 117   Double_t resolFull=resol->GetFullEvResol();
 118   Double_t resolFull2=AliEventPlaneResolution::GetFullEvResol(hResolSubAB,1);
 119   printf("Resolution on full event = %.4f %.4f\n",resolFull,resolFull2);
 120
 121   TCanvas* cEP=new TCanvas("cEP","EvPlaneResol");
 122   hResolSubAB->Draw();
 123   TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
 124   tres->SetNDC();
 125   tres->Draw();
 126   cEP->Update();
 127
 128   TCanvas** c1=new TCanvas*[nFinalPtBins];
 129   TCanvas** c2=new TCanvas*[nFinalPtBins];
 130   TF1** fSig=new TF1*[nFinalPtBins];
 131   TF1** fv2=new TF1*[nFinalPtBins];
 132   TH1F** hAverCos2Phi=new TH1F*[nFinalPtBins];
 133   TH1F** hFractionSig=new TH1F*[nFinalPtBins];
 134   TH1F** hFractionBkg=new TH1F*[nFinalPtBins];
 135
 136   TLegend* leg0=0x0;
 137   TLegend* leg1=0x0;
 138   TLegendEntry* ent;
 139
 140   TH1F* hMass=0x0;
 141   TH1F* hCos2PhiBkgLo=0x0;
 142   TH1F* hCos2PhiBkgHi=0x0;
 143   TH1F* hCos2PhiBkgLoScal=0x0;
 144   TH1F* hCos2PhiBkgHiScal=0x0;
 145   TH1F* hCos2PhiBkgAver=0x0;
 146   TH1F* hCos2PhiSigReg=0x0;
 147   TH1F* hCos2PhiSig=0x0;
 148
 149   Double_t v2M1[nFinalPtBins],errv2M1[nFinalPtBins];
 150   Double_t v2M2[nFinalPtBins],errv2M2[nFinalPtBins];
 151
 152   for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 153     printf("**************** BIN %d ******************\n",iFinalPtBin);
 154     printf("\n--------- Method 1: Side Bands ----------\n");
 155     hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
 156
 157     Double_t sigmaRangeForSig=3.;
 158     Double_t sigmaRangeForBkg=4.5;
 159
 160     gStyle->SetPalette(1);
 161     gStyle->SetOptTitle(0);
 162     c1[iFinalPtBin]=new TCanvas(Form("cMeth1PtBin%d",iFinalPtBin),Form("Meth1PtBin%d",iFinalPtBin));
 163     c1[iFinalPtBin]->Divide(2,2);
 164     c1[iFinalPtBin]->cd(1);
 165     hMassDphi[iFinalPtBin]->Draw("colz");
 166     c1[iFinalPtBin]->cd(2);
 167
 168     Int_t nMassBins=hMass->GetNbinsX();
 169     Int_t hMinBin=3;
 170     Int_t hMaxBin=nMassBins-2;
 171     Double_t hmin=hMass->GetBinLowEdge(hMinBin);
 172     Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
 173     Int_t factor4refl=0;
 174     Float_t massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 175
 176     AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
 177     fitter->SetReflectionSigmaFactor(factor4refl);
 178     fitter->SetInitialGaussianMean(massD);
 179     Bool_t out=fitter->MassFitter(0);
 180     if(!out) continue;
 181     fitter->DrawHere(gPad);
 182     Double_t sigfitter,esigfitter;
 183     fitter->Signal(sigmaRangeForSig, sigfitter,esigfitter);
 184     Double_t mass=fitter->GetMean();
 185     Double_t sigma=fitter->GetSigma();
 186     TF1* fB1=fitter->GetBackgroundFullRangeFunc();
 187     TF1* fB2=fitter->GetBackgroundRecalcFunc();
 188     TF1* fSB=fitter->GetMassFunc();
 189     Double_t minMassSig=mass-sigmaRangeForSig*sigma;
 190     Double_t maxMassSig=mass+sigmaRangeForSig*sigma;
 191     Int_t minBinSig=hMass->FindBin(minMassSig);
 192     Int_t maxBinSig=hMass->FindBin(maxMassSig);
 193     Double_t minMassSigBin=hMass->GetBinLowEdge(minBinSig);
 194     Double_t maxMassSigBin=hMass->GetBinLowEdge(maxBinSig)+hMass->GetBinWidth(maxBinSig);
 195     printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
 196     Double_t maxMassBkgLow=mass-sigmaRangeForBkg*sigma;
 197     Int_t minBinBkgLow=hMinBin;
 198     Int_t maxBinBkgLow=hMass->FindBin(maxMassBkgLow);
 199     Double_t minMassBkgLowBin=hmin;
 200     Double_t maxMassBkgLowBin=hMass->GetBinLowEdge(maxBinBkgLow)+hMass->GetBinWidth(maxBinBkgLow);
 201     Double_t minMassBkgHi=mass+sigmaRangeForBkg*sigma;
 202     Int_t minBinBkgHi=hMass->FindBin(minMassBkgHi);
 203     Int_t maxBinBkgHi=hMaxBin;
 204     Double_t minMassBkgHiBin=hMass->GetBinLowEdge(minBinBkgHi);
 205     Double_t maxMassBkgHiBin=hmax;
 206     printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
 207     Double_t bkgSig=fB2->Integral(minMassSigBin,maxMassSigBin);
 208     Double_t bkgLow=fB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
 209     Double_t bkgHi=fB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
 210     printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
 211     TBox* bleft=new TBox(minMassBkgLowBin,0.,maxMassBkgLowBin,hMass->GetMaximum());
 212     bleft->SetFillColor(kRed+1);
 213     bleft->SetFillStyle(3002);
 214     bleft->Draw();
 215     TBox* bright=new TBox(minMassBkgHiBin,0.,maxMassBkgHiBin,hMass->GetMaximum());
 216     bright->SetFillColor(kBlue+1);
 217     bright->SetFillStyle(3002);
 218     bright->Draw();
 219     TBox* bsig=new TBox(minMassSigBin,0.,maxMassSigBin,hMass->GetMaximum()*2);
 220     bsig->SetFillColor(1);
 221     bsig->SetFillStyle(3002);
 222     bsig->Draw();
 223
 224     hCos2PhiBkgLo=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionX(Form("hCos2PhiBkgLoBin%d",iFinalPtBin),minBinBkgLow,maxBinBkgLow);
 225     hCos2PhiBkgHi=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionX(Form("hCos2PhiBkgHiBin%d",iFinalPtBin),minBinBkgHi,maxBinBkgHi);
 226     hCos2PhiSigReg=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionX(Form("hCos2PhiBkgSigBin%d",iFinalPtBin),minBinSig,maxBinSig);
 227
 228     hCos2PhiBkgLo->Rebin(4);
 229     hCos2PhiBkgHi->Rebin(4);
 230     hCos2PhiSigReg->Rebin(4);
 231     hCos2PhiSigReg->SetLineWidth(2);
 232     hCos2PhiBkgLo->SetLineWidth(2);
 233     hCos2PhiBkgHi->SetLineWidth(2);
 234     hCos2PhiBkgLo->SetLineColor(kRed+1);
 235     hCos2PhiBkgHi->SetLineColor(kBlue+1);
 236
 237     hCos2PhiBkgLoScal=(TH1F*)hCos2PhiBkgLo->Clone(Form("hCos2PhiBkgLoScalBin%d",iFinalPtBin));
 238     hCos2PhiBkgLoScal->Scale(bkgSig/bkgLow);
 239     hCos2PhiBkgHiScal=(TH1F*)hCos2PhiBkgHi->Clone(Form("hCos2PhiBkgHiScalBin%d",iFinalPtBin));
 240     hCos2PhiBkgHiScal->Scale(bkgSig/bkgHi);
 241     hCos2PhiBkgLoScal->SetLineWidth(2);
 242     hCos2PhiBkgHiScal->SetLineWidth(2);
 243     hCos2PhiBkgLoScal->SetLineColor(kRed+1);
 244     hCos2PhiBkgHiScal->SetLineColor(kBlue+1);
 245     hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
 246     hCos2PhiBkgAver->Add(hCos2PhiBkgHiScal);
 247     hCos2PhiBkgAver->Scale(0.5);
 248     hCos2PhiBkgAver->SetLineWidth(2);
 249     hCos2PhiBkgAver->SetLineColor(kGreen+1);
 250     hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone(Form("hCos2PhiSigBin%d",iFinalPtBin));
 251     hCos2PhiSig->Add(hCos2PhiBkgAver,-1.);
 252
 253     c1[iFinalPtBin]->cd(3);
 254     hCos2PhiSigReg->Draw();
 255     hCos2PhiBkgLoScal->Draw("same");
 256     hCos2PhiBkgHiScal->Draw("same");
 257     hCos2PhiBkgAver->Draw("same");
 258     if(iFinalPtBin==0){
 259       leg0=new TLegend(0.3,0.6,0.75,0.89);
 260       leg0->SetFillColor(0);
 261       ent=leg0->AddEntry(hCos2PhiSigReg,"Signal region","L");
 262       ent->SetTextColor(hCos2PhiSigReg->GetLineColor());
 263       ent=leg0->AddEntry(hCos2PhiBkgLoScal,"Left side band","L");
 264       ent->SetTextColor(hCos2PhiBkgLoScal->GetLineColor());
 265       ent=leg0->AddEntry(hCos2PhiBkgHiScal,"Right side band","L");
 266       ent->SetTextColor(hCos2PhiBkgHiScal->GetLineColor());
 267       ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","L");
 268       ent->SetTextColor(hCos2PhiBkgAver->GetLineColor());
 269     }
 270     leg0->Draw();
 271     c1[iFinalPtBin]->cd(4);
 272     hCos2PhiSig->Draw("EP");
 273     TPaveText* t0= new TPaveText(0.15,0.70,0.45,0.89,"NDC");
 274     t0->SetFillColor(0);
 275     t0->AddText(Form("v2=%.3f+-%.3f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError()));
 276     t0->Draw();
 277
 278     printf("Signal from mass fitter = %f  Signal from subracted histo= %f\n",
 279            sigfitter,hCos2PhiSig->Integral());
 280     Double_t v2obsM1=hCos2PhiSig->GetMean();
 281     Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
 282     printf("v2obs = %f +- %f\n",v2obsM1,errv2obsM1);
 283     v2M1[iFinalPtBin]=v2obsM1/resolFull;
 284     errv2M1[iFinalPtBin]=errv2obsM1/resolFull;
 285     printf("v2 = %f +- %f\n",v2M1[iFinalPtBin],errv2M1[iFinalPtBin]);
 286
 287     printf("\n--------- Method 2: S/S+B ----------\n");
 288     Int_t npars=fSB->GetNpar();
 289     Double_t sigmaSB=fSB->GetParameter(npars-1);
 290     Double_t massSB=fSB->GetParameter(npars-2);
 291     Double_t integr=fSB->GetParameter(npars-3);
 292     Double_t sOverAll=(fSB->Eval(massSB)-fB2->Eval(massSB))/fSB->Eval(massSB);
 293     printf("mass=%f  S+B=%f   bkg=%f S/(S+B)=%f\n",massSB,fSB->Eval(massSB),fB2->Eval(massSB),sOverAll);
 294     printf("Number of parameters: %d. Signal params: %f %f %f\n",npars,massSB,sigmaSB,integr);
 295
 296     Int_t nbinsmass=hMassDphi[iFinalPtBin]->GetNbinsY();
 297     Double_t minmass=hMassDphi[iFinalPtBin]->GetYaxis()->GetXmin();
 298     Double_t maxmass=hMassDphi[iFinalPtBin]->GetYaxis()->GetXmax();
 299
 300     fSig[iFinalPtBin]=new TF1(Form("fSig%d",iFinalPtBin),"[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./[2]/[2])",minmass,maxmass);
 301     fSig[iFinalPtBin]->SetParameters(integr,massSB,sigmaSB);
 302
 303     hAverCos2Phi[iFinalPtBin]=new TH1F(Form("hAverCos2PhiBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 304     hFractionSig[iFinalPtBin]=new TH1F(Form("hFractionSigBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 305     hFractionBkg[iFinalPtBin]=new TH1F(Form("hFractionBkgBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
 306
 307     for(Int_t iBin=1; iBin<=hMassDphi[iFinalPtBin]->GetNbinsY(); iBin++){
 308       TH1F* htemp=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionX("htemp",iBin,iBin);
 309       hAverCos2Phi[iFinalPtBin]->SetBinContent(iBin,htemp->GetMean());
 310       hAverCos2Phi[iFinalPtBin]->SetBinError(iBin,htemp->GetMeanError());
 311       Double_t sig=fSig[iFinalPtBin]->Eval(hFractionSig[iFinalPtBin]->GetBinCenter(iBin));
 312       Double_t bkg=fB2->Eval(hFractionSig[iFinalPtBin]->GetBinCenter(iBin));
 313       if(bkg<1 && sig<1){
 314         hFractionSig[iFinalPtBin]->SetBinContent(iBin,0.);
 315         hFractionSig[iFinalPtBin]->SetBinError(iBin,0.);
 316         hFractionBkg[iFinalPtBin]->SetBinContent(iBin,1.);
 317         hFractionBkg[iFinalPtBin]->SetBinError(iBin,0.);
 318       }else{
 319         Double_t fracs=sig/(sig+bkg);
 320         Double_t fracb=bkg/(sig+bkg);
 321         Double_t efracs=0.;//TMath::Sqrt(fracs*(1.-fracs)/(sig+bkg));
 322         Double_t efracb=0.;//TMath::Sqrt(fracb*(1.-fracb)/(sig+bkg));
 323
 324         hFractionSig[iFinalPtBin]->SetBinContent(iBin,fracs);
 325         hFractionSig[iFinalPtBin]->SetBinError(iBin,efracs);
 326         hFractionBkg[iFinalPtBin]->SetBinContent(iBin,fracb);
 327         hFractionBkg[iFinalPtBin]->SetBinError(iBin,efracb);
 328       }
 329       delete htemp;
 330     }
 331
 332     fv2[iFinalPtBin]=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMass,minmass,maxmass,5);
 333     fv2[iFinalPtBin]->SetParameter(0,sOverAll);
 334     fv2[iFinalPtBin]->SetParameter(1,massSB);
 335     fv2[iFinalPtBin]->SetParameter(2,sigmaSB);
 336     fv2[iFinalPtBin]->SetParameter(3,0.2);
 337     fv2[iFinalPtBin]->SetParameter(4,0.2);
 338     fv2[iFinalPtBin]->FixParameter(0,sOverAll);
 339     fv2[iFinalPtBin]->FixParameter(1,massSB);
 340     fv2[iFinalPtBin]->FixParameter(2,sigmaSB);
 341
 342     hAverCos2Phi[iFinalPtBin]->Rebin(2);
 343     hAverCos2Phi[iFinalPtBin]->Scale(0.5);
 344
 345     c2[iFinalPtBin]=new TCanvas(Form("cMeth2Bin%d",iFinalPtBin),Form("cMeth2Bin%d",iFinalPtBin));
 346     c2[iFinalPtBin]->Divide(2,2);
 347     c2[iFinalPtBin]->cd(1);
 348     hMassDphi[iFinalPtBin]->Draw("colz");
 349     c2[iFinalPtBin]->cd(2);
 350     hMass->Rebin(2);
 351     hMass->SetMinimum(0.);
 352     hMass->SetMarkerStyle(20);
 353     hMass->Draw("E");
 354     fSB->Draw("same");
 355     fSig[iFinalPtBin]->Draw("same");
 356     fB2->Draw("same");
 357     c2[iFinalPtBin]->cd(3);
 358     hFractionSig[iFinalPtBin]->SetMaximum(1.2);
 359     hFractionSig[iFinalPtBin]->Draw();
 360     hFractionSig[iFinalPtBin]->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 361     hFractionSig[iFinalPtBin]->GetYaxis()->SetTitle("Fraction");
 362     hFractionBkg[iFinalPtBin]->SetLineColor(2);
 363     hFractionBkg[iFinalPtBin]->Draw("same");
 364     if(iFinalPtBin==0){
 365       leg1=new TLegend(0.15,0.15,0.35,0.35);
 366       leg1->SetFillColor(0);
 367       ent=leg1->AddEntry(hFractionSig[iFinalPtBin],"S/(S+B)","L");
 368       ent->SetTextColor(hFractionSig[iFinalPtBin]->GetLineColor());
 369       ent=leg1->AddEntry(hFractionBkg[iFinalPtBin],"B/(S+B)","L");
 370       ent->SetTextColor(hFractionBkg[iFinalPtBin]->GetLineColor());
 371     }
 372     leg1->Draw();
 373     c2[iFinalPtBin]->cd(4);
 374     hAverCos2Phi[iFinalPtBin]->Fit(fv2[iFinalPtBin]);
 375     hAverCos2Phi[iFinalPtBin]->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 376     hAverCos2Phi[iFinalPtBin]->GetYaxis()->SetTitle("v_2^{obs}");
 377     TPaveText* t1= new TPaveText(0.55,0.70,0.89,0.89,"NDC");
 378     t1->SetFillColor(0);
 379     t1->AddText(Form("v2sig=%.3f+-%.3f\n",fv2[iFinalPtBin]->GetParameter(3),fv2[iFinalPtBin]->GetParError(3)));
 380     t1->AddText(Form("v2bkg=%.3f+-%.3f\n",fv2[iFinalPtBin]->GetParameter(4),fv2[iFinalPtBin]->GetParError(4)));
 381     t1->Draw();
 382
 383     Double_t v2obsM2=fv2[iFinalPtBin]->GetParameter(3);
 384     Double_t errv2obsM2=fv2[iFinalPtBin]->GetParError(3);
 385     printf("v2obs = %f +- %f\n",v2obsM2,errv2obsM2);
 386     v2M2[iFinalPtBin]=v2obsM2/resolFull;
 387     errv2M2[iFinalPtBin]=errv2obsM2/resolFull;
 388     printf("v2 = %f +- %f\n",v2M2[iFinalPtBin],errv2M2[iFinalPtBin]);
 389   }
 390
 391   printf("\n--------- Summary ------------\n");
 392
 393   TH1F* hv2m1=new TH1F("hv2m1","",nFinalPtBins,ptlims);
 394   TH1F* hv2m2=new TH1F("hv2m2","",nFinalPtBins,ptlims);
 395    for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
 396     printf("PtBin %d   v2method1 = %f +- %f   v2method2 = %f +-%f\n",iFinalPtBin,
 397            v2M1[iFinalPtBin],errv2M1[iFinalPtBin],
 398            v2M2[iFinalPtBin],errv2M2[iFinalPtBin]
 399            );
 400     hv2m1->SetBinContent(iFinalPtBin+1,v2M1[iFinalPtBin]);
 401     hv2m1->SetBinError(iFinalPtBin+1,errv2M1[iFinalPtBin]);
 402     hv2m2->SetBinContent(iFinalPtBin+1,v2M2[iFinalPtBin]);
 403     hv2m2->SetBinError(iFinalPtBin+1,errv2M2[iFinalPtBin]);
 404   }
 405
 406    Double_t maxy=TMath::Max(hv2m2->GetMaximum(),hv2m1->GetMaximum())+0.1;
 407    Double_t miny=TMath::Min(hv2m2->GetMinimum(),hv2m1->GetMinimum())-0.1;
 408    TH2F* hempty=new TH2F("hempty","",10,0.,hv2m1->GetXaxis()->GetXmax()+2.,10,miny,maxy);
 409    hempty->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 410    hempty->GetYaxis()->SetTitle("v_{2}");
 411
 412    TCanvas* cv2=new TCanvas("cv2","v2");
 413    hempty->Draw();
 414    hv2m1->SetMarkerStyle(26);
 415    hv2m1->Draw("same");
 416    hv2m2->SetLineColor(2);
 417    hv2m2->SetMarkerColor(2);
 418    hv2m2->SetMarkerStyle(20);
 419    hv2m2->Draw("same");
 420    TLegend* leg2=new TLegend(0.5,0.7,0.89,0.89);
 421    leg2->SetFillStyle(0);
 422    ent=leg2->AddEntry(hv2m1,"Side Band subtraction","P");
 423    ent->SetTextColor(hv2m1->GetMarkerColor());
 424    ent=leg2->AddEntry(hv2m2,"Fit to v2 vs. mass","P");
 425    ent->SetTextColor(hv2m2->GetMarkerColor());
 426    leg2->Draw();
 427    cv2->Update();
 428 }