]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/charmFlow/Extractv2from2Dhistos.C
Updates in the macro for v2 calculation from 2D (mass vs. cos2phi) histos
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / charmFlow / Extractv2from2Dhistos.C
CommitLineData
30c6e405 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
a8f6c03f 27
28Double_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
41void Extractv2from2Dhistos(){
42
30c6e405 43 gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
44
45 TFile* fil=new TFile("AnalysisResults_ptbins.root");
a8f6c03f 46 TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
30c6e405 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
a8f6c03f 55 Int_t minCent=30;
56 Int_t maxCent=50;
30c6e405 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
a8f6c03f 82 for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){
30c6e405 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 }
a8f6c03f 106 }
107
30c6e405 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;
a8f6c03f 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;
a8f6c03f 148
30c6e405 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));
a8f6c03f 323
30c6e405 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;
a8f6c03f 330 }
a8f6c03f 331
30c6e405 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}");
a8f6c03f 411
30c6e405 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();
a8f6c03f 428}