| 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> |
| 07515723 |
19 | #include <TGraphErrors.h> |
| 20 | #include <TGraphAsymmErrors.h> |
| 21 | |
| 30c6e405 |
22 | #include "AliAODRecoDecayHF.h" |
| 23 | #include "AliRDHFCuts.h" |
| 24 | #include "AliRDHFCutsDplustoKpipi.h" |
| 98d04a62 |
25 | #include "AliRDHFCutsDStartoKpipi.h" |
| 30c6e405 |
26 | #include "AliRDHFCutsD0toKpi.h" |
| 27 | #include "AliHFMassFitter.h" |
| 0784bdc7 |
28 | #include "AliEventPlaneResolutionHandler.h" |
| 5d540e24 |
29 | #include "AliVertexingHFUtils.h" |
| 30c6e405 |
30 | #endif |
| a8f6c03f |
31 | |
| a8f6c03f |
32 | |
| 07515723 |
33 | // Common variables: to be configured by the user |
| 0784bdc7 |
34 | TString filename="AnalysisResults_train63.root"; |
| 30d8161c |
35 | TString suffix="v2Dplus3050Cut4upcutPIDTPC"; |
| 36 | TString partname="Dplus"; |
| 07515723 |
37 | Int_t minCent=30; |
| 38 | Int_t maxCent=50; |
| 0784bdc7 |
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; |
| 6760e3ca |
46 | |
| 30d8161c |
47 | const Int_t nFinalPtBins=4; |
| 48 | Double_t ptlims[nFinalPtBins+1]={3.,4.,6.,8.,12.}; |
| 07515723 |
49 | Double_t sigmaRangeForSig=2.5; |
| 50 | Double_t sigmaRangeForBkg=4.5; |
| 30d8161c |
51 | Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2,2}; |
| 6760e3ca |
52 | Bool_t useConstantvsBkgVsMass=kFALSE; |
| 30d8161c |
53 | Int_t rebinHistov2Mass[nFinalPtBins]={2,2,2,2}; |
| 07515723 |
54 | Int_t factor4refl=0; |
| 30d8161c |
55 | Int_t minPtBin[nFinalPtBins]={-1,-1,-1,-1}; |
| 56 | Int_t maxPtBin[nFinalPtBins]={-1,-1,-1,-1}; |
| 98d04a62 |
57 | Bool_t saveAllCanvas=kFALSE; |
| 30d8161c |
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 | }; |
| d94072a7 |
73 | |
| 6760e3ca |
74 | // systematic errors for 2-5, 5-8 and 8-12 (no run-by-run weights) |
| 75 | /* |
| 07515723 |
76 | Double_t systErrMeth1[nFinalPtBins]={ |
| d94072a7 |
77 | (0.308-0.169)/2., |
| 78 | (0.14-0.1)/2., |
| 79 | (0.04+0.02)/2. |
| 07515723 |
80 | }; |
| 81 | Double_t systErrMeth2[nFinalPtBins]={ |
| d94072a7 |
82 | (0.305-0.252)/2., |
| 83 | (0.129-0.020)/2., |
| 84 | (0.101+0.06)/2. |
| 07515723 |
85 | }; |
| 6760e3ca |
86 | */ |
| 07515723 |
87 | |
| 6760e3ca |
88 | // systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights) |
| 30d8161c |
89 | /* |
| 6760e3ca |
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. |
| 0784bdc7 |
99 | ; |
| 30d8161c |
100 | */ |
| 101 | |
| 6760e3ca |
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 | */ |
| 07515723 |
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); |
| 98d04a62 |
131 | Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3); |
| 07515723 |
132 | Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad); |
| 6760e3ca |
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 | |
| a8f6c03f |
136 | |
| 98d04a62 |
137 | //______________________________________________________________________________ |
| a8f6c03f |
138 | void Extractv2from2Dhistos(){ |
| 30d8161c |
139 | // main function: computes v2 with side band and v2(M) methods |
| a8f6c03f |
140 | |
| 30c6e405 |
141 | gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C"); |
| 142 | |
| 30d8161c |
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); |
| 08580a2e |
157 | if(!binOK){ |
| 158 | printf("ERROR: mismatch in pt binning\n"); |
| 159 | return; |
| 160 | } |
| 30d8161c |
161 | TList *lst =(TList*)dir->Get(listname.Data()); |
| 6760e3ca |
162 | if(!lst){ |
| 30d8161c |
163 | printf("list %s not found in file, please check list name\n",listname.Data()); |
| 6760e3ca |
164 | return; |
| 165 | } |
| 30d8161c |
166 | |
| 0784bdc7 |
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); |
| 30c6e405 |
181 | |
| 30c6e405 |
182 | TH2F** hMassDphi=new TH2F*[nFinalPtBins]; |
| 183 | for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){ |
| 184 | hMassDphi[iFinalPtBin]=0x0; |
| 185 | } |
| 07515723 |
186 | LoadMassHistos(lst, hMassDphi); |
| 30c6e405 |
187 | |
| 188 | TCanvas** c1=new TCanvas*[nFinalPtBins]; |
| 189 | TCanvas** c2=new TCanvas*[nFinalPtBins]; |
| 07515723 |
190 | |
| 30c6e405 |
191 | Double_t v2M1[nFinalPtBins],errv2M1[nFinalPtBins]; |
| 192 | Double_t v2M2[nFinalPtBins],errv2M2[nFinalPtBins]; |
| 193 | |
| 07515723 |
194 | gStyle->SetPalette(1); |
| 195 | gStyle->SetOptTitle(0); |
| 196 | |
| 30c6e405 |
197 | for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){ |
| 198 | printf("**************** BIN %d ******************\n",iFinalPtBin); |
| 199 | printf("\n--------- Method 1: Side Bands ----------\n"); |
| 07515723 |
200 | TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY(); |
| 30c6e405 |
201 | |
| 30c6e405 |
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); |
| 07515723 |
207 | Bool_t out=FitMassSpectrum(hMass,(TPad*)gPad); |
| 30c6e405 |
208 | if(!out) continue; |
| 30c6e405 |
209 | |
| 6760e3ca |
210 | TH1F* hCos2PhiSig=DoSideBands(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],c1[iFinalPtBin]); |
| 30c6e405 |
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"); |
| 04140e26 |
219 | c2[iFinalPtBin]=new TCanvas(Form("cMeth2PtBin%d",iFinalPtBin),Form("Meth2PtBin%d",iFinalPtBin)); |
| 6760e3ca |
220 | TF1* fv2=DoFitv2vsMass(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],c2[iFinalPtBin],0,useConstantvsBkgVsMass); |
| 30c6e405 |
221 | |
| 07515723 |
222 | Double_t v2obsM2=fv2->GetParameter(3); |
| 223 | Double_t errv2obsM2=fv2->GetParError(3); |
| 30c6e405 |
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]); |
| 98d04a62 |
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 | } |
| 30c6e405 |
234 | } |
| 235 | |
| 236 | printf("\n--------- Summary ------------\n"); |
| 237 | |
| 8460a05f |
238 | TH1F* hv2m1=new TH1F("hv2m1","Side Band subtraction",nFinalPtBins,ptlims); |
| 239 | TH1F* hv2m2=new TH1F("hv2m2","Fit to v2 vs. mass",nFinalPtBins,ptlims); |
| 30c6e405 |
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 | |
| 07515723 |
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]; |
| d94072a7 |
256 | Double_t syste1=TMath::Sqrt(systErrMeth1[iFinalPtBin]*systErrMeth1[iFinalPtBin]+systRes1*systRes1); |
| 07515723 |
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 | |
| 30c6e405 |
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}"); |
| a8f6c03f |
268 | |
| 30c6e405 |
269 | TCanvas* cv2=new TCanvas("cv2","v2"); |
| 270 | hempty->Draw(); |
| 271 | hv2m1->SetMarkerStyle(26); |
| 07515723 |
272 | hSystErr1->SetFillColor(kGray); |
| 273 | hSystErr1->SetFillStyle(3002); |
| 274 | hSystErr1->Draw("e2same"); |
| 275 | |
| 276 | hv2m2->SetLineColor(kRed+1); |
| 277 | hv2m2->SetMarkerColor(kRed+1); |
| 30c6e405 |
278 | hv2m2->SetMarkerStyle(20); |
| 07515723 |
279 | hSystErr2->SetFillColor(kRed-9); |
| 280 | hSystErr2->SetFillStyle(3005); |
| 281 | hSystErr2->Draw("e2same"); |
| 282 | |
| 283 | hv2m1->Draw("same"); |
| 30c6e405 |
284 | hv2m2->Draw("same"); |
| 07515723 |
285 | |
| 30c6e405 |
286 | TLegend* leg2=new TLegend(0.5,0.7,0.89,0.89); |
| 287 | leg2->SetFillStyle(0); |
| 07515723 |
288 | TLegendEntry* ent=leg2->AddEntry(hv2m1,"Side Band subtraction","P"); |
| 30c6e405 |
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(); |
| 98d04a62 |
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())); |
| 8460a05f |
296 | |
| 30d8161c |
297 | TFile* outfil=new TFile(Form("v2Output2DMeth_%s_%d_%d_%s.root",partname.Data(),minCent,maxCent,suffix.Data()),"recreate"); |
| 8460a05f |
298 | outfil->cd(); |
| 299 | hv2m1->Write(); |
| 300 | hv2m2->Write(); |
| d94072a7 |
301 | hSystErr1->Write(); |
| 302 | hSystErr2->Write(); |
| 8460a05f |
303 | outfil->Close(); |
| a8f6c03f |
304 | } |
| 07515723 |
305 | |
| 306 | |
| 98d04a62 |
307 | //______________________________________________________________________________ |
| 07515723 |
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 | |
| 98d04a62 |
321 | //______________________________________________________________________________ |
| d94072a7 |
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 | |
| 98d04a62 |
336 | //______________________________________________________________________________ |
| 07515723 |
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); |
| 98d04a62 |
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 | } |
| 07515723 |
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 | |
| 98d04a62 |
376 | //______________________________________________________________________________ |
| 07515723 |
377 | TH1F* DoSideBands(Int_t iFinalPtBin, |
| 378 | TH2F* hMassDphi, |
| 379 | TH1F* hMass, |
| 6760e3ca |
380 | Int_t rebin, |
| d94072a7 |
381 | TCanvas* c1, |
| 382 | Int_t optBkg){ |
| 07515723 |
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); |
| 04140e26 |
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 | |
| 6760e3ca |
437 | hCos2PhiBkgLo->Rebin(rebin); |
| 438 | hCos2PhiBkgHi->Rebin(rebin); |
| 439 | hCos2PhiSigReg->Rebin(rebin); |
| 07515723 |
440 | hCos2PhiSigReg->SetLineWidth(2); |
| 441 | hCos2PhiBkgLo->SetLineWidth(2); |
| 442 | hCos2PhiBkgHi->SetLineWidth(2); |
| 443 | hCos2PhiBkgLo->SetLineColor(kRed+1); |
| 444 | hCos2PhiBkgHi->SetLineColor(kBlue+1); |
| 04140e26 |
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(); |
| 07515723 |
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); |
| 04140e26 |
462 | hCos2PhiBkgLoScal->SetLineStyle(7); |
| 463 | hCos2PhiBkgHiScal->SetLineStyle(2); |
| 07515723 |
464 | hCos2PhiBkgLoScal->SetLineColor(kRed+1); |
| 465 | hCos2PhiBkgHiScal->SetLineColor(kBlue+1); |
| 04140e26 |
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 | |
| d94072a7 |
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 | } |
| 04140e26 |
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 | |
| 07515723 |
487 | hCos2PhiBkgAver->SetLineWidth(2); |
| 488 | hCos2PhiBkgAver->SetLineColor(kGreen+1); |
| 04140e26 |
489 | hCos2PhiBkgAver->SetFillColor(kGreen+1); |
| 07515723 |
490 | TH1F* hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone(Form("hCos2PhiSigBin%d",iFinalPtBin)); |
| 491 | hCos2PhiSig->Add(hCos2PhiBkgAver,-1.); |
| 04140e26 |
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 | ); |
| 07515723 |
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(); |
| 04140e26 |
503 | hCos2PhiBkgAver->Draw("same"); |
| 07515723 |
504 | hCos2PhiBkgLoScal->Draw("same"); |
| 505 | hCos2PhiBkgHiScal->Draw("same"); |
| 07515723 |
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()); |
| 04140e26 |
513 | ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","F"); |
| 07515723 |
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 | |
| 98d04a62 |
540 | //______________________________________________________________________________ |
| d94072a7 |
541 | TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors, Bool_t useConst){ |
| 07515723 |
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); |
| d94072a7 |
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 | } |
| 07515723 |
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); |
| 07515723 |
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 | |
| d94072a7 |
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 | } |
| 07515723 |
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 | |
| 30d8161c |
608 | printf("LOOP %d\n",rebin); |
| d94072a7 |
609 | if((hAverCos2Phi->GetNbinsX()%rebin)==0){ |
| 610 | hAverCos2Phi->Rebin(rebin); |
| 611 | hAverCos2Phi->Scale(1./(Double_t)rebin); |
| 07515723 |
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); |
| d94072a7 |
642 | fv2->DrawCopy("same"); |
| 07515723 |
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))); |
| d94072a7 |
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 | } |
| 07515723 |
653 | t1->Draw(); |
| 654 | c2->Update(); |
| 655 | }else{ |
| 656 | hAverCos2Phi->Fit(fv2); |
| 657 | } |
| 658 | return fv2; |
| 659 | } |
| 660 | |
| 661 | |
| 98d04a62 |
662 | //______________________________________________________________________________ |
| 07515723 |
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; |
| 98d04a62 |
666 | TH1F* hResolSubBC=0x0; |
| 667 | TH1F* hResolSubAC=0x0; |
| 07515723 |
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; |
| 30d8161c |
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){ |
| 98d04a62 |
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()); |
| 30d8161c |
681 | cout<<hResolSubABsing->GetName()<<endl; |
| 98d04a62 |
682 | TH1F* hResolSubBCsing=0x0; |
| 683 | TH1F* hResolSubACsing=0x0; |
| 0784bdc7 |
684 | if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub|| |
| 685 | evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){ |
| 98d04a62 |
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; |
| 30d8161c |
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); |
| 07515723 |
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; |
| 30d8161c |
705 | if(iHisC==minCentTimesTen){ |
| 07515723 |
706 | hResolSubAB=(TH1F*)hResolSubABsing->Clone("hResolSubAB"); |
| 0784bdc7 |
707 | if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub|| |
| 708 | evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){ |
| 98d04a62 |
709 | hResolSubBC=(TH1F*)hResolSubBCsing->Clone("hResolSubAB"); |
| 710 | hResolSubAC=(TH1F*)hResolSubACsing->Clone("hResolSubAB"); |
| 711 | } |
| 07515723 |
712 | }else{ |
| 713 | hResolSubAB->Add(hResolSubABsing); |
| 0784bdc7 |
714 | if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub|| |
| 715 | evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){ |
| 98d04a62 |
716 | hResolSubBC->Add(hResolSubBCsing); |
| 717 | hResolSubAC->Add(hResolSubACsing); |
| 718 | } |
| 07515723 |
719 | } |
| 98d04a62 |
720 | printf("Adding histogram %s\n",hResolSubAB->GetName()); |
| 07515723 |
721 | } |
| 722 | rcflow=xmin; |
| 723 | rcfhigh=xmax; |
| 724 | |
| 5d540e24 |
725 | |
| 98d04a62 |
726 | Double_t error; |
| 727 | Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubAB,hResolSubBC,hResolSubAC); |
| 07515723 |
728 | |
| 729 | grInteg->SetPoint(0,40.,resolFull); |
| 730 | grInteg->SetPointEXlow(0,10); |
| 731 | grInteg->SetPointEXhigh(0,10.); |
| 98d04a62 |
732 | grInteg->SetPointEYlow(0,error); |
| 733 | grInteg->SetPointEYhigh(0,error); |
| 07515723 |
734 | |
| 07515723 |
735 | TCanvas* cEP=new TCanvas("cEP","EvPlaneResol"); |
| 736 | cEP->Divide(1,2); |
| 737 | cEP->cd(1); |
| 738 | hResolSubAB->Draw(); |
| 0784bdc7 |
739 | if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub|| |
| 740 | evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){ |
| 98d04a62 |
741 | hResolSubBC->SetLineColor(2); |
| 742 | hResolSubBC->Draw("same"); |
| 743 | hResolSubAC->SetLineColor(4); |
| 744 | hResolSubAC->Draw("same"); |
| 745 | } |
| 07515723 |
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 | |
| 98d04a62 |
763 | //______________________________________________________________________________ |
| 764 | Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){ |
| 765 | Double_t resolFull; |
| 0784bdc7 |
766 | if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub || |
| 767 | evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){ |
| 98d04a62 |
768 | resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1); |
| 769 | error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1)); |
| 0784bdc7 |
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 | } |
| 98d04a62 |
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]); |
| 0784bdc7 |
790 | if(evPlane==AliEventPlaneResolutionHandler::kVZEROC || |
| 791 | evPlane==AliEventPlaneResolutionHandler::kVZERO){ |
| 98d04a62 |
792 | resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]); |
| 793 | error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]); |
| 794 | } |
| 0784bdc7 |
795 | else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){ |
| 98d04a62 |
796 | resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]); |
| 797 | error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]); |
| 798 | } |
| 0784bdc7 |
799 | else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta || |
| 800 | evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){ |
| 98d04a62 |
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 | } |
| 07515723 |
807 | |
| 98d04a62 |
808 | //______________________________________________________________________________ |
| 07515723 |
809 | void LoadMassHistos(TList* lst, TH2F** hMassDphi){ |
| 810 | |
| 30d8161c |
811 | Int_t minCentTimesTen=minCent*10; |
| 812 | Int_t maxCentTimesTen=maxCent*10; |
| 813 | for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){ |
| 07515723 |
814 | for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){ |
| 815 | for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){ |
| 30d8161c |
816 | TString hisname=Form("hMc2deltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25); |
| 07515723 |
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 | |
| 98d04a62 |
829 | //______________________________________________________________________________ |
| 07515723 |
830 | Bool_t DefinePtBins(TDirectoryFile* df){ |
| 98d04a62 |
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 | |
| 07515723 |
837 | Int_t nPtBinsCuts=d0cuts->GetNPtBins(); |
| 838 | printf("Number of pt bins for cut object = %d\n",nPtBinsCuts); |
| 839 | Float_t *ptlimsCuts=d0cuts->GetPtBinLimits(); |
| 08580a2e |
840 | for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]); |
| 07515723 |
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; |
| 08580a2e |
845 | if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1; |
| 07515723 |
846 | } |
| 847 | } |
| 08580a2e |
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; |
| 07515723 |
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 | } |
| 08580a2e |
855 | |
| 07515723 |
856 | return kTRUE; |
| 857 | } |
| 858 | |
| 859 | |
| 98d04a62 |
860 | //______________________________________________________________________________ |
| 07515723 |
861 | void SystForSideBands(){ |
| 30d8161c |
862 | // Compute systematic ucnertainty for the side band subtraction method |
| 07515723 |
863 | gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C"); |
| 864 | |
| 30d8161c |
865 | TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data()); |
| 866 | TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data()); |
| d94072a7 |
867 | |
| 30d8161c |
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()); |
| 6760e3ca |
884 | if(!lst){ |
| 30d8161c |
885 | printf("list %s not found in file, please check list name\n",listname.Data()); |
| 6760e3ca |
886 | return; |
| 887 | } |
| 07515723 |
888 | |
| 0784bdc7 |
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); |
| 07515723 |
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]; |
| d94072a7 |
908 | TGraphErrors** gSystWhichSide=new TGraphErrors*[nFinalPtBins]; |
| 07515723 |
909 | |
| 910 | Double_t min1[nFinalPtBins],max1[nFinalPtBins]; |
| 6760e3ca |
911 | Double_t min123[nFinalPtBins],max123[nFinalPtBins]; |
| 07515723 |
912 | Double_t min2[nFinalPtBins],max2[nFinalPtBins]; |
| 913 | Double_t min3[nFinalPtBins],max3[nFinalPtBins]; |
| d94072a7 |
914 | Double_t min4[nFinalPtBins],max4[nFinalPtBins]; |
| 07515723 |
915 | |
| 916 | Double_t sigmaRangeForSigOrig=sigmaRangeForSig; |
| 917 | Double_t sigmaRangeForBkgOrig=sigmaRangeForBkg; |
| 07515723 |
918 | |
| 919 | |
| 920 | for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){ |
| 921 | printf("**************** BIN %d ******************\n",iFinalPtBin); |
| 922 | |
| 6760e3ca |
923 | Int_t rebinHistoSideBandsOrig=rebinHistoSideBands[iFinalPtBin]; |
| 924 | |
| 07515723 |
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.; |
| 6760e3ca |
935 | min123[iFinalPtBin]=99999.; |
| 936 | max123[iFinalPtBin]=-99999.; |
| 07515723 |
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; |
| 6760e3ca |
941 | rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig; |
| 942 | TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0); |
| 07515723 |
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; |
| 6760e3ca |
952 | if(sigmaRangeForSig>=2. && sigmaRangeForSig<=3){ |
| 953 | if(v2M1>max123[iFinalPtBin]) max123[iFinalPtBin]=v2M1; |
| 954 | if(v2M1<min123[iFinalPtBin]) min123[iFinalPtBin]=v2M1; |
| 955 | } |
| 07515723 |
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; |
| 6760e3ca |
966 | rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig; |
| 967 | TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0); |
| 07515723 |
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; |
| 6760e3ca |
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); |
| 07515723 |
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; |
| 6760e3ca |
996 | gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistoSideBands[iFinalPtBin],v2M1); |
| 07515723 |
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 | } |
| d94072a7 |
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; |
| 6760e3ca |
1012 | rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig; |
| 1013 | TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0,iStep); |
| d94072a7 |
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 | } |
| 6760e3ca |
1025 | rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig; |
| 07515723 |
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]); |
| 6760e3ca |
1031 | printf(" (limited to 2-3 sigma) = %f %f\n",min123[iFinalPtBin],max123[iFinalPtBin]); |
| 07515723 |
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]); |
| d94072a7 |
1034 | printf("Range of values for whichside = %f %f\n",min4[iFinalPtBin],max4[iFinalPtBin]); |
| 6760e3ca |
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); |
| 07515723 |
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 | } |
| d94072a7 |
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 | } |
| 07515723 |
1079 | } |
| 1080 | |
| 98d04a62 |
1081 | //______________________________________________________________________________ |
| 07515723 |
1082 | void SystForFitv2Mass(){ |
| 30d8161c |
1083 | // Compute systematic ucnertainty for the v2(M) method |
| 07515723 |
1084 | gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C"); |
| 1085 | |
| 30d8161c |
1086 | TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data()); |
| 1087 | TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data()); |
| d94072a7 |
1088 | |
| 30d8161c |
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()); |
| 6760e3ca |
1105 | if(!lst){ |
| 30d8161c |
1106 | printf("list %s not found in file, please check list name\n",listname.Data()); |
| 6760e3ca |
1107 | return; |
| 1108 | } |
| 07515723 |
1109 | |
| 0784bdc7 |
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 | |
| 07515723 |
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 | |
| d94072a7 |
1126 | TGraphErrors** gSystParamErr=new TGraphErrors*[nFinalPtBins]; |
| 07515723 |
1127 | TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins]; |
| d94072a7 |
1128 | TGraphErrors** gSystLinConst=new TGraphErrors*[nFinalPtBins]; |
| 07515723 |
1129 | |
| 1130 | Double_t min1[nFinalPtBins],max1[nFinalPtBins]; |
| d94072a7 |
1131 | Double_t min2[nFinalPtBins],max2[nFinalPtBins]; |
| 1132 | Double_t min3[nFinalPtBins],max3[nFinalPtBins]; |
| 07515723 |
1133 | |
| 07515723 |
1134 | |
| 1135 | for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){ |
| 1136 | printf("**************** BIN %d ******************\n",iFinalPtBin); |
| 1137 | |
| d94072a7 |
1138 | Int_t rebinHistov2MassOrig=rebinHistov2Mass[iFinalPtBin]; |
| 07515723 |
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; |
| d94072a7 |
1151 | rebinHistov2Mass[iFinalPtBin]=iStep+1; |
| 1152 | if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistov2Mass[iFinalPtBin]!=0) continue; |
| 6760e3ca |
1153 | TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,useConstantvsBkgVsMass); |
| 07515723 |
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; |
| d94072a7 |
1159 | gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistov2Mass[iFinalPtBin],v2M2); |
| 07515723 |
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 | } |
| d94072a7 |
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; |
| 6760e3ca |
1174 | TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,iStep,useConstantvsBkgVsMass); |
| d94072a7 |
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 | |
| 07515723 |
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]); |
| d94072a7 |
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]); |
| 6760e3ca |
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); |
| 07515723 |
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 | } |
| d94072a7 |
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 | } |
| 07515723 |
1250 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame