+++ /dev/null
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <TInterpreter.h>
-#include <TString.h>
-#include <TObjString.h>
-#include <TObjArray.h>
-#include <TMath.h>
-#include <TFile.h>
-#include <TCanvas.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TH1D.h>
-#include <TF1.h>
-#include <TStyle.h>
-#include <TLegend.h>
-#include <TLegendEntry.h>
-#include <TLatex.h>
-#include <TPaveText.h>
-#include <TDatabasePDG.h>
-#include <TGraphErrors.h>
-#include <TGraphAsymmErrors.h>
-
-#include "AliAODRecoDecayHF.h"
-#include "AliRDHFCuts.h"
-#include "AliRDHFCutsDplustoKpipi.h"
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliRDHFCutsD0toKpi.h"
-#include "AliHFMassFitter.h"
-#include "AliEventPlaneResolutionHandler.h"
-#include "AliVertexingHFUtils.h"
-#endif
-
-
-// Common variables: to be configured by the user
-TString filename="AnalysisResults_train63.root";
-TString suffix="v2Dplus3050Cut4upcutPIDTPC";
-TString partname="Dplus";
-Int_t minCent=30;
-Int_t maxCent=50;
-//evPlane flag from AliEventPlaneResolutionHandler:
-//kTPCFullEta, kTPCPosEta,kVZERO,kVZEROA,kVZEROC
-Int_t evPlane=AliEventPlaneResolutionHandler::kTPCPosEta;
-//resolution flag fromAliEventPlaneResolutionHandler:
-//kTwoRandSub,kTwoChargeSub,kTwoEtaSub,kThreeSub,kThreeSubTPCGap
-Int_t evPlaneRes=AliEventPlaneResolutionHandler::kTwoEtaSub;
-Bool_t useNcollWeight=kFALSE;
-
-const Int_t nFinalPtBins=4;
-Double_t ptlims[nFinalPtBins+1]={3.,4.,6.,8.,12.};
-Double_t sigmaRangeForSig=2.5;
-Double_t sigmaRangeForBkg=4.5;
-Int_t rebinHistoSideBands[nFinalPtBins]={2,2,2,2};
-Bool_t useConstantvsBkgVsMass=kFALSE;
-Int_t rebinHistov2Mass[nFinalPtBins]={2,2,2,2};
-Int_t factor4refl=0;
-Int_t minPtBin[nFinalPtBins]={-1,-1,-1,-1};
-Int_t maxPtBin[nFinalPtBins]={-1,-1,-1,-1};
-Bool_t saveAllCanvas=kFALSE;
-
-// systematic errors for D+, 30-50 TPC eventplane
-Double_t systErrMeth1[nFinalPtBins]={
- (0.247359-0.152503)/2.,
- (0.337073-0.290978)/2.,
- (0.286888-0.222045)/2.,
- (-0.078888+0.178259)/2.
-};
-
-Double_t systErrMeth2[nFinalPtBins]={
- (0.151019-0.066624)/2.,
- (0.272252-0.232044)/2.,
- (0.246019-0.143002)/2.,
- (-0.026077+0.292956)/2.
-};
-
-// systematic errors for 2-5, 5-8 and 8-12 (no run-by-run weights)
-/*
-Double_t systErrMeth1[nFinalPtBins]={
- (0.308-0.169)/2.,
- (0.14-0.1)/2.,
- (0.04+0.02)/2.
-};
-Double_t systErrMeth2[nFinalPtBins]={
- (0.305-0.252)/2.,
- (0.129-0.020)/2.,
- (0.101+0.06)/2.
-};
-*/
-
-// systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights)
-/*
-Double_t systErrMeth1[nFinalPtBins]={
- (0.23-0.10)/2.,
- (0.152-0.078)/2.,
- (0.161-0.097)/2.
-};
-Double_t systErrMeth2[nFinalPtBins]={
- (0.164-0.097)/2.,
- (0.110-0.012)/2.,
- (0.131-0.036)/2.
-;
-*/
-
-
-
-// systematic errors for 2-5, 5-8 and 8-12 (93 runs with run-by-run weights, RAA cuts)
-/*
-Double_t systErrMeth1[nFinalPtBins]={
- (0.265-0.122)/2.,
- (0.165-0.117)/2.,
- (0.238-0.169)/2.
-};
-Double_t systErrMeth2[nFinalPtBins]={
- (0.174-0.135)/2.,
- (0.18-0.11)/2.,
- (0.311-0.28)/2.
-};
-*/
-
-// output of mass fitter
-Int_t hMinBin;
-Int_t hMaxBin;
-Double_t signalFromFit,esignalFromFit;
-Double_t massFromFit,sigmaFromFit;
-TF1* fB1=0x0;
-TF1* fB2=0x0;
-TF1* fSB=0x0;
-
-void LoadMassHistos(TList* lst, TH2F** hMassDphi);
-Bool_t DefinePtBins(TDirectoryFile* df);
-Double_t v2vsMass(Double_t *x, Double_t *par);
-Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh);
-Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);
-Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad);
-TH1F* DoSideBands(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c1, Int_t optBkg=0);
-TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors=0, Bool_t useConst=kTRUE);
-
-
-//______________________________________________________________________________
-void Extractv2from2Dhistos(){
- // main function: computes v2 with side band and v2(M) methods
-
- gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
-
- TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
- TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
-
- TFile *f = TFile::Open(filename.Data());
- if(!f){
- printf("file %s not found, please check file name\n",filename.Data());
- return;
- }
- TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
- if(!dir){
- printf("Directory %s not found, please check dir name\n",dirname.Data());
- return;
- }
- Bool_t binOK=DefinePtBins(dir);
- if(!binOK){
- printf("ERROR: mismatch in pt binning\n");
- return;
- }
- TList *lst =(TList*)dir->Get(listname.Data());
- if(!lst){
- printf("list %s not found in file, please check list name\n",listname.Data());
- return;
- }
-
- // Double_t rcfmin,rcfmax;
- // Double_t resolFull=GetEPResolution(lst,rcfmin,rcfmax);
- // Double_t resolSyst=(rcfmax-rcfmin)/2./resolFull;
- // printf("Relative Systematic Error on RCF=%f\n",resolSyst);
- AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
- epres->SetEventPlane(evPlane);
- epres->SetResolutionOption(evPlaneRes);
- epres->SetUseNcollWeights(useNcollWeight);
- Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
- Double_t rcfmin=epres->GetEventPlaneResolution(minCent,minCent+2.5);
- Double_t rcfmax=epres->GetEventPlaneResolution(maxCent-2.5,maxCent);
- Double_t resolSyst=TMath::Abs(rcfmax-rcfmin)/2./resolFull;
- delete epres;
- printf("Event plane resolution %f\n",resolFull);
-
- TH2F** hMassDphi=new TH2F*[nFinalPtBins];
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- hMassDphi[iFinalPtBin]=0x0;
- }
- LoadMassHistos(lst, hMassDphi);
-
- TCanvas** c1=new TCanvas*[nFinalPtBins];
- TCanvas** c2=new TCanvas*[nFinalPtBins];
-
- Double_t v2M1[nFinalPtBins],errv2M1[nFinalPtBins];
- Double_t v2M2[nFinalPtBins],errv2M2[nFinalPtBins];
-
- gStyle->SetPalette(1);
- gStyle->SetOptTitle(0);
-
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("**************** BIN %d ******************\n",iFinalPtBin);
- printf("\n--------- Method 1: Side Bands ----------\n");
- TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
-
- c1[iFinalPtBin]=new TCanvas(Form("cMeth1PtBin%d",iFinalPtBin),Form("Meth1PtBin%d",iFinalPtBin));
- c1[iFinalPtBin]->Divide(2,2);
- c1[iFinalPtBin]->cd(1);
- hMassDphi[iFinalPtBin]->Draw("colz");
- c1[iFinalPtBin]->cd(2);
- Bool_t out=FitMassSpectrum(hMass,(TPad*)gPad);
- if(!out) continue;
-
- TH1F* hCos2PhiSig=DoSideBands(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],c1[iFinalPtBin]);
- Double_t v2obsM1=hCos2PhiSig->GetMean();
- Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
- printf("v2obs = %f +- %f\n",v2obsM1,errv2obsM1);
- v2M1[iFinalPtBin]=v2obsM1/resolFull;
- errv2M1[iFinalPtBin]=errv2obsM1/resolFull;
- printf("v2 = %f +- %f\n",v2M1[iFinalPtBin],errv2M1[iFinalPtBin]);
-
- printf("\n--------- Method 2: S/S+B ----------\n");
- c2[iFinalPtBin]=new TCanvas(Form("cMeth2PtBin%d",iFinalPtBin),Form("Meth2PtBin%d",iFinalPtBin));
- TF1* fv2=DoFitv2vsMass(iFinalPtBin,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],c2[iFinalPtBin],0,useConstantvsBkgVsMass);
-
- Double_t v2obsM2=fv2->GetParameter(3);
- Double_t errv2obsM2=fv2->GetParError(3);
- printf("v2obs = %f +- %f\n",v2obsM2,errv2obsM2);
- v2M2[iFinalPtBin]=v2obsM2/resolFull;
- errv2M2[iFinalPtBin]=errv2obsM2/resolFull;
- printf("v2 = %f +- %f\n",v2M2[iFinalPtBin],errv2M2[iFinalPtBin]);
- if(saveAllCanvas){
- c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.root",iFinalPtBin));
- c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.root",iFinalPtBin));
- c1[iFinalPtBin]->SaveAs(Form("cMeth1PtBin%d.eps",iFinalPtBin));
- c2[iFinalPtBin]->SaveAs(Form("cMeth2PtBin%d.eps",iFinalPtBin));
- }
- }
-
- printf("\n--------- Summary ------------\n");
-
- TH1F* hv2m1=new TH1F("hv2m1","Side Band subtraction",nFinalPtBins,ptlims);
- TH1F* hv2m2=new TH1F("hv2m2","Fit to v2 vs. mass",nFinalPtBins,ptlims);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("PtBin %d v2method1 = %f +- %f v2method2 = %f +-%f\n",iFinalPtBin,
- v2M1[iFinalPtBin],errv2M1[iFinalPtBin],
- v2M2[iFinalPtBin],errv2M2[iFinalPtBin]
- );
- hv2m1->SetBinContent(iFinalPtBin+1,v2M1[iFinalPtBin]);
- hv2m1->SetBinError(iFinalPtBin+1,errv2M1[iFinalPtBin]);
- hv2m2->SetBinContent(iFinalPtBin+1,v2M2[iFinalPtBin]);
- hv2m2->SetBinError(iFinalPtBin+1,errv2M2[iFinalPtBin]);
- }
-
- TH1F* hSystErr1=(TH1F*)hv2m1->Clone("hSystErr1");
- TH1F* hSystErr2=(TH1F*)hv2m2->Clone("hSystErr2");
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- Double_t systRes1=resolSyst*v2M1[iFinalPtBin];
- Double_t systRes2=resolSyst*v2M2[iFinalPtBin];
- Double_t syste1=TMath::Sqrt(systErrMeth1[iFinalPtBin]*systErrMeth1[iFinalPtBin]+systRes1*systRes1);
- Double_t syste2=TMath::Sqrt(systErrMeth2[iFinalPtBin]*systErrMeth2[iFinalPtBin]+systRes2*systRes2);
- hSystErr1->SetBinError(iFinalPtBin+1,syste1);
- hSystErr2->SetBinError(iFinalPtBin+1,syste2);
- }
-
-
- Double_t maxy=TMath::Max(hv2m2->GetMaximum(),hv2m1->GetMaximum())+0.1;
- Double_t miny=TMath::Min(hv2m2->GetMinimum(),hv2m1->GetMinimum())-0.1;
- TH2F* hempty=new TH2F("hempty","",10,0.,hv2m1->GetXaxis()->GetXmax()+2.,10,miny,maxy);
- hempty->GetXaxis()->SetTitle("p_{t} (GeV/c)");
- hempty->GetYaxis()->SetTitle("v_{2}");
-
- TCanvas* cv2=new TCanvas("cv2","v2");
- hempty->Draw();
- hv2m1->SetMarkerStyle(26);
- hSystErr1->SetFillColor(kGray);
- hSystErr1->SetFillStyle(3002);
- hSystErr1->Draw("e2same");
-
- hv2m2->SetLineColor(kRed+1);
- hv2m2->SetMarkerColor(kRed+1);
- hv2m2->SetMarkerStyle(20);
- hSystErr2->SetFillColor(kRed-9);
- hSystErr2->SetFillStyle(3005);
- hSystErr2->Draw("e2same");
-
- hv2m1->Draw("same");
- hv2m2->Draw("same");
-
- TLegend* leg2=new TLegend(0.5,0.7,0.89,0.89);
- leg2->SetFillStyle(0);
- TLegendEntry* ent=leg2->AddEntry(hv2m1,"Side Band subtraction","P");
- ent->SetTextColor(hv2m1->GetMarkerColor());
- ent=leg2->AddEntry(hv2m2,"Fit to v2 vs. mass","P");
- ent->SetTextColor(hv2m2->GetMarkerColor());
- leg2->Draw();
- cv2->Update();
- cv2->SaveAs(Form("v2%s-2Dmethods-%s.png",partname.Data(),suffix.Data()));
- cv2->SaveAs(Form("v2%s-2Dmethods-%s.eps",partname.Data(),suffix.Data()));
-
- TFile* outfil=new TFile(Form("v2Output2DMeth_%s_%d_%d_%s.root",partname.Data(),minCent,maxCent,suffix.Data()),"recreate");
- outfil->cd();
- hv2m1->Write();
- hv2m2->Write();
- hSystErr1->Write();
- hSystErr2->Write();
- outfil->Close();
-}
-
-
-//______________________________________________________________________________
-Double_t v2vsMass(Double_t *x, Double_t *par){
- // Fit function for signal+background
- // par[0] = S/B at the mass peak
- // par[1] = mass
- // par[2] = sigma
- // par[3] = v2sig
- // par[4] = v2back
-
- Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
- Double_t fracbkg=1-fracsig;
- return par[3]*fracsig+par[4]*fracbkg;
-}
-
-//______________________________________________________________________________
-Double_t v2vsMassLin(Double_t *x, Double_t *par){
- // Fit function for signal+background
- // par[0] = S/B at the mass peak
- // par[1] = mass
- // par[2] = sigma
- // par[3] = v2sig
- // par[4] = v2back constant
- // par[5] = v2back slope
-
- Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
- Double_t fracbkg=1-fracsig;
- return par[3]*fracsig+(par[4]+par[5]*x[0])*fracbkg;
-}
-
-//______________________________________________________________________________
-Bool_t FitMassSpectrum(TH1F* hMass, TPad* pad){
-
- Int_t nMassBins=hMass->GetNbinsX();
- hMinBin=3;
- hMaxBin=nMassBins-2;
- Double_t hmin=hMass->GetBinLowEdge(hMinBin);
- Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
- Double_t massD;
- if(partname.Contains("Dzero")) {
- massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
- }else if(partname.Contains("Dplus")){
- massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
- }else if(partname.Contains("Dstar")) {
- massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass());
- }else{
- printf("Wrong particle name\n");
- massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
- }
-
- AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
- fitter->SetReflectionSigmaFactor(factor4refl);
- fitter->SetInitialGaussianMean(massD);
- Bool_t out=fitter->MassFitter(0);
- if(!out) return kFALSE;
- fitter->Signal(sigmaRangeForSig, signalFromFit,esignalFromFit);
- massFromFit=fitter->GetMean();
- sigmaFromFit=fitter->GetSigma();
- fB1=fitter->GetBackgroundFullRangeFunc();
- fB2=fitter->GetBackgroundRecalcFunc();
- fSB=fitter->GetMassFunc();
- if(!fB1) return kFALSE;
- if(!fB2) return kFALSE;
- if(!fSB) return kFALSE;
- if(pad){
- fitter->DrawHere(gPad,3.,0.);
- }
- return kTRUE;
-}
-
-//______________________________________________________________________________
-TH1F* DoSideBands(Int_t iFinalPtBin,
- TH2F* hMassDphi,
- TH1F* hMass,
- Int_t rebin,
- TCanvas* c1,
- Int_t optBkg){
-
- // Build histo with cos(2*deltaphi) distribution for signal
-
- Double_t hmin=hMass->GetBinLowEdge(hMinBin);
- Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
-
- Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
- Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
- Int_t minBinSig=hMass->FindBin(minMassSig);
- Int_t maxBinSig=hMass->FindBin(maxMassSig);
- Double_t minMassSigBin=hMass->GetBinLowEdge(minBinSig);
- Double_t maxMassSigBin=hMass->GetBinLowEdge(maxBinSig)+hMass->GetBinWidth(maxBinSig);
- printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
- Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
- Int_t minBinBkgLow=hMinBin;
- Int_t maxBinBkgLow=hMass->FindBin(maxMassBkgLow);
- Double_t minMassBkgLowBin=hmin;
- Double_t maxMassBkgLowBin=hMass->GetBinLowEdge(maxBinBkgLow)+hMass->GetBinWidth(maxBinBkgLow);
- Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
- Int_t minBinBkgHi=hMass->FindBin(minMassBkgHi);
- Int_t maxBinBkgHi=hMaxBin;
- Double_t minMassBkgHiBin=hMass->GetBinLowEdge(minBinBkgHi);
- Double_t maxMassBkgHiBin=hmax;
- printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
- Double_t bkgSig=fB2->Integral(minMassSigBin,maxMassSigBin);
- Double_t bkgLow=fB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
- Double_t bkgHi=fB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
- printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
- if(c1){
- TBox* bleft=new TBox(minMassBkgLowBin,0.,maxMassBkgLowBin,hMass->GetMaximum());
- bleft->SetFillColor(kRed+1);
- bleft->SetFillStyle(3002);
- bleft->Draw();
- TBox* bright=new TBox(minMassBkgHiBin,0.,maxMassBkgHiBin,hMass->GetMaximum());
- bright->SetFillColor(kBlue+1);
- bright->SetFillStyle(3002);
- bright->Draw();
- TBox* bsig=new TBox(minMassSigBin,0.,maxMassSigBin,hMass->GetMaximum()*2);
- bsig->SetFillColor(1);
- bsig->SetFillStyle(3002);
- bsig->Draw();
- }
- TH1F* hCos2PhiBkgLo=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgLoBin%d",iFinalPtBin),minBinBkgLow,maxBinBkgLow);
- TH1F* hCos2PhiBkgHi=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgHiBin%d",iFinalPtBin),minBinBkgHi,maxBinBkgHi);
- TH1F* hCos2PhiSigReg=(TH1F*)hMassDphi->ProjectionX(Form("hCos2PhiBkgSigBin%d",iFinalPtBin),minBinSig,maxBinSig);
-
- printf("Before Rebin11: BkgLo = %f+-%f BkgHi=%f+-%f SigReg=%f+-%f\n",hCos2PhiBkgLo->GetBinContent(11),hCos2PhiBkgLo->GetBinError(11),
- hCos2PhiBkgHi->GetBinContent(11),hCos2PhiBkgHi->GetBinError(11),
- hCos2PhiSigReg->GetBinContent(11),hCos2PhiSigReg->GetBinError(11));
- printf("Before Rebin12: BkgLo = %f+-%f BkgHi=%f+-%f SigReg=%f+-%f\n",hCos2PhiBkgLo->GetBinContent(12),hCos2PhiBkgLo->GetBinError(12),
- hCos2PhiBkgHi->GetBinContent(12),hCos2PhiBkgHi->GetBinError(12),
- hCos2PhiSigReg->GetBinContent(12),hCos2PhiSigReg->GetBinError(12));
- Double_t binc=hCos2PhiBkgLo->GetBinCenter(11);
-
- hCos2PhiBkgLo->Rebin(rebin);
- hCos2PhiBkgHi->Rebin(rebin);
- hCos2PhiSigReg->Rebin(rebin);
- hCos2PhiSigReg->SetLineWidth(2);
- hCos2PhiBkgLo->SetLineWidth(2);
- hCos2PhiBkgHi->SetLineWidth(2);
- hCos2PhiBkgLo->SetLineColor(kRed+1);
- hCos2PhiBkgHi->SetLineColor(kBlue+1);
-
- Int_t theBinR=hCos2PhiBkgLo->FindBin(binc);
- printf("After Rebin %d: BkgLo = %f+-%f BkgHi=%f+-%f SigReg=%f+-%f\n",theBinR,
- hCos2PhiBkgLo->GetBinContent(theBinR),hCos2PhiBkgLo->GetBinError(theBinR),
- hCos2PhiBkgHi->GetBinContent(theBinR),hCos2PhiBkgHi->GetBinError(theBinR),
- hCos2PhiSigReg->GetBinContent(theBinR),hCos2PhiSigReg->GetBinError(theBinR));
-
- printf("Background Scaling: Lo->Sig= %f Hi->Sig=%f\n",bkgSig/bkgLow,bkgSig/bkgHi);
-
- hCos2PhiBkgLo->Sumw2();
- hCos2PhiBkgHi->Sumw2();
- TH1F* hCos2PhiBkgLoScal=(TH1F*)hCos2PhiBkgLo->Clone(Form("hCos2PhiBkgLoScalBin%d",iFinalPtBin));
- hCos2PhiBkgLoScal->Scale(bkgSig/bkgLow);
- TH1F* hCos2PhiBkgHiScal=(TH1F*)hCos2PhiBkgHi->Clone(Form("hCos2PhiBkgHiScalBin%d",iFinalPtBin));
- hCos2PhiBkgHiScal->Scale(bkgSig/bkgHi);
- hCos2PhiBkgLoScal->SetLineWidth(2);
- hCos2PhiBkgHiScal->SetLineWidth(2);
- hCos2PhiBkgLoScal->SetLineStyle(7);
- hCos2PhiBkgHiScal->SetLineStyle(2);
- hCos2PhiBkgLoScal->SetLineColor(kRed+1);
- hCos2PhiBkgHiScal->SetLineColor(kBlue+1);
- printf("After Scaling %d: BkgLo = %f+-%f BkgHi=%f+-%f \n",theBinR,
- hCos2PhiBkgLoScal->GetBinContent(theBinR),hCos2PhiBkgLoScal->GetBinError(theBinR),
- hCos2PhiBkgHiScal->GetBinContent(theBinR),hCos2PhiBkgHiScal->GetBinError(theBinR));
-
-
- TH1F* hCos2PhiBkgAver=0x0;
- if(optBkg==0){
- hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
- hCos2PhiBkgAver->Add(hCos2PhiBkgHiScal);
- hCos2PhiBkgAver->Scale(0.5);
- }else if(optBkg==-1){
- hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
- }else{
- hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgHiScal->Clone(Form("hCos2PhiBkgAverBin%d",iFinalPtBin));
- }
- printf("After Average %d: BkgLo = %f+-%f BkgHi=%f+-%f BkgAve=%f+-%f \n",theBinR,
- hCos2PhiBkgLoScal->GetBinContent(theBinR),hCos2PhiBkgLoScal->GetBinError(theBinR),
- hCos2PhiBkgHiScal->GetBinContent(theBinR),hCos2PhiBkgHiScal->GetBinError(theBinR),
- hCos2PhiBkgAver->GetBinContent(theBinR),hCos2PhiBkgAver->GetBinError(theBinR)
- );
-
- hCos2PhiBkgAver->SetLineWidth(2);
- hCos2PhiBkgAver->SetLineColor(kGreen+1);
- hCos2PhiBkgAver->SetFillColor(kGreen+1);
- TH1F* hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone(Form("hCos2PhiSigBin%d",iFinalPtBin));
- hCos2PhiSig->Add(hCos2PhiBkgAver,-1.);
- printf("After Subtraction %d: All = %f+-%f BkgAve=%f+-%f Sig=%f+-%f \n",theBinR,
- hCos2PhiSigReg->GetBinContent(theBinR),hCos2PhiSigReg->GetBinError(theBinR),
- hCos2PhiBkgAver->GetBinContent(theBinR),hCos2PhiBkgAver->GetBinError(theBinR),
- hCos2PhiSig->GetBinContent(theBinR),hCos2PhiSig->GetBinError(theBinR)
- );
-
- TLegend* leg0=new TLegend(0.3,0.6,0.75,0.89);
- TPaveText* t0= new TPaveText(0.15,0.70,0.45,0.89,"NDC");
- if(c1){
- c1->cd(3);
- hCos2PhiSigReg->Draw();
- hCos2PhiBkgAver->Draw("same");
- hCos2PhiBkgLoScal->Draw("same");
- hCos2PhiBkgHiScal->Draw("same");
- leg0->SetFillColor(0);
- TLegendEntry* ent=leg0->AddEntry(hCos2PhiSigReg,"Signal region","L");
- ent->SetTextColor(hCos2PhiSigReg->GetLineColor());
- ent=leg0->AddEntry(hCos2PhiBkgLoScal,"Left side band","L");
- ent->SetTextColor(hCos2PhiBkgLoScal->GetLineColor());
- ent=leg0->AddEntry(hCos2PhiBkgHiScal,"Right side band","L");
- ent->SetTextColor(hCos2PhiBkgHiScal->GetLineColor());
- ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","F");
- ent->SetTextColor(hCos2PhiBkgAver->GetLineColor());
- leg0->Draw();
- c1->cd(4);
- hCos2PhiSig->Draw("EP");
- t0->SetFillColor(0);
- t0->AddText(Form("v2=%.3f+-%.3f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError()));
- t0->Draw();
-
- c1->Update();
- }
- if(!c1){
- delete leg0;
- delete t0;
- delete hCos2PhiBkgLo;
- delete hCos2PhiBkgHi;
- delete hCos2PhiSigReg;
- delete hCos2PhiBkgLoScal;
- delete hCos2PhiBkgHiScal;
- delete hCos2PhiBkgAver;
- }
- printf("Signal from mass fitter = %f Signal from subracted histo= %f\n",
- signalFromFit,hCos2PhiSig->Integral());
- return hCos2PhiSig;
-}
-
-
-//______________________________________________________________________________
-TF1* DoFitv2vsMass(Int_t iFinalPtBin, TH2F* hMassDphi, TH1F* hMass, Int_t rebin, TCanvas* c2, Int_t optErrors, Bool_t useConst){
-
- Int_t npars=fSB->GetNpar();
- Double_t sigmaSB=fSB->GetParameter(npars-1);
- Double_t massSB=fSB->GetParameter(npars-2);
- Double_t integr=fSB->GetParameter(npars-3);
- Double_t sOverAll=(fSB->Eval(massSB)-fB2->Eval(massSB))/fSB->Eval(massSB);
- printf("mass=%f S+B=%f bkg=%f S/(S+B)=%f\n",massSB,fSB->Eval(massSB),fB2->Eval(massSB),sOverAll);
- printf("Number of parameters: %d. Signal params: %f %f %f\n",npars,massSB,sigmaSB,integr);
- if(optErrors==1){
- sOverAll+=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
- }
- else if(optErrors==-1){
- sOverAll-=fSB->GetParError(npars-3)/fSB->GetParameter(npars-3);
- }
-
- Int_t nbinsmass=hMassDphi->GetNbinsY();
- Double_t minmass=hMassDphi->GetYaxis()->GetXmin();
- Double_t maxmass=hMassDphi->GetYaxis()->GetXmax();
-
- 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);
- fSig->SetParameters(integr,massSB,sigmaSB);
-
- TH1F* hAverCos2Phi=new TH1F(Form("hAverCos2PhiBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
- TH1F* hFractionSig=new TH1F(Form("hFractionSigBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
- TH1F* hFractionBkg=new TH1F(Form("hFractionBkgBin%d",iFinalPtBin),"",nbinsmass,minmass,maxmass);
- for(Int_t iBin=1; iBin<=hMassDphi->GetNbinsY(); iBin++){
- TH1F* htemp=(TH1F*)hMassDphi->ProjectionX("htemp",iBin,iBin);
- hAverCos2Phi->SetBinContent(iBin,htemp->GetMean());
- hAverCos2Phi->SetBinError(iBin,htemp->GetMeanError());
- Double_t sig=fSig->Eval(hFractionSig->GetBinCenter(iBin));
- Double_t bkg=fB2->Eval(hFractionSig->GetBinCenter(iBin));
- if(bkg<1 && sig<1){
- hFractionSig->SetBinContent(iBin,0.);
- hFractionSig->SetBinError(iBin,0.);
- hFractionBkg->SetBinContent(iBin,1.);
- hFractionBkg->SetBinError(iBin,0.);
- }else{
- Double_t fracs=sig/(sig+bkg);
- Double_t fracb=bkg/(sig+bkg);
- Double_t efracs=0.;//TMath::Sqrt(fracs*(1.-fracs)/(sig+bkg));
- Double_t efracb=0.;//TMath::Sqrt(fracb*(1.-fracb)/(sig+bkg));
-
- hFractionSig->SetBinContent(iBin,fracs);
- hFractionSig->SetBinError(iBin,efracs);
- hFractionBkg->SetBinContent(iBin,fracb);
- hFractionBkg->SetBinError(iBin,efracb);
- }
- delete htemp;
- }
-
- TF1* fv2=0x0;
- if(useConst){
- fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMass,minmass,maxmass,5);
- }else{
- fv2=new TF1(Form("fv2Bin%d",iFinalPtBin),v2vsMassLin,minmass,maxmass,6);
- fv2->SetParameter(5,0.);
- }
- fv2->SetParameter(0,sOverAll);
- fv2->SetParameter(1,massSB);
- fv2->SetParameter(2,sigmaSB);
- fv2->SetParameter(3,0.2);
- fv2->SetParameter(4,0.2);
- fv2->FixParameter(0,sOverAll);
- fv2->FixParameter(1,massSB);
- fv2->FixParameter(2,sigmaSB);
-
- printf("LOOP %d\n",rebin);
- if((hAverCos2Phi->GetNbinsX()%rebin)==0){
- hAverCos2Phi->Rebin(rebin);
- hAverCos2Phi->Scale(1./(Double_t)rebin);
- }
-
- if(c2){
- c2->Divide(2,2);
- c2->cd(1);
- hMassDphi->Draw("colz");
- c2->cd(2);
- hMass->Rebin(2);
- hMass->SetMinimum(0.);
- hMass->SetMarkerStyle(20);
- hMass->Draw("E");
- fSB->Draw("same");
- fSig->Draw("same");
- fB2->Draw("same");
- c2->cd(3);
- hFractionSig->SetMaximum(1.2);
- hFractionSig->Draw();
- hFractionSig->GetXaxis()->SetTitle("Mass (GeV/c^2)");
- hFractionSig->GetYaxis()->SetTitle("Fraction");
- hFractionBkg->SetLineColor(2);
- hFractionBkg->Draw("same");
- TLegend* leg1=new TLegend(0.15,0.15,0.35,0.35);
- leg1->SetFillColor(0);
- TLegendEntry* ent=leg1->AddEntry(hFractionSig,"S/(S+B)","L");
- ent->SetTextColor(hFractionSig->GetLineColor());
- ent=leg1->AddEntry(hFractionBkg,"B/(S+B)","L");
- ent->SetTextColor(hFractionBkg->GetLineColor());
- leg1->Draw();
- c2->cd(4);
- hAverCos2Phi->Fit(fv2);
- fv2->DrawCopy("same");
- hAverCos2Phi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
- hAverCos2Phi->GetYaxis()->SetTitle("v_2^{obs}");
- TPaveText* t1= new TPaveText(0.55,0.70,0.89,0.89,"NDC");
- t1->SetFillColor(0);
- t1->AddText(Form("v2sig=%.3f+-%.3f\n",fv2->GetParameter(3),fv2->GetParError(3)));
- if(useConst){
- t1->AddText(Form("v2bkg=%.3f+-%.3f\n",fv2->GetParameter(4),fv2->GetParError(4)));
- }else{
- t1->AddText(Form("v2bkg=(%.3f+-%.3f) + (%.3g+-%.3g)*mass\n",fv2->GetParameter(4),fv2->GetParError(4),fv2->GetParameter(5),fv2->GetParError(5)));
- }
- t1->Draw();
- c2->Update();
- }else{
- hAverCos2Phi->Fit(fv2);
- }
- return fv2;
-}
-
-
-//______________________________________________________________________________
-Double_t GetEPResolution(TList* lst, Double_t &rcflow, Double_t &rcfhigh){
- // Event plane resolution syst err (from wide centrality bin
- TH1F* hResolSubAB=0x0;
- TH1F* hResolSubBC=0x0;
- TH1F* hResolSubAC=0x0;
- Double_t xmin=1.;
- Double_t xmax=-1.;
- TGraphAsymmErrors* grSingle=new TGraphAsymmErrors(0);
- TGraphAsymmErrors* grInteg=new TGraphAsymmErrors(0);
- Int_t iPt=0;
- Int_t minCentTimesTen=minCent*10;
- Int_t maxCentTimesTen=maxCent*10;
- Double_t binHalfWid=25./20.;
- for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){
- TString hisnameAB=Form("hEvPlaneResocentr%d_%d",iHisC,iHisC+25);
- TString hisnameBC=Form("hEvPlaneReso2centr%d_%d",iHisC,iHisC+25);
- TString hisnameAC=Form("hEvPlaneReso3centr%d_%d",iHisC,iHisC+25);
- TH1F* hResolSubABsing=(TH1F*)lst->FindObject(hisnameAB.Data());
- cout<<hResolSubABsing->GetName()<<endl;
- TH1F* hResolSubBCsing=0x0;
- TH1F* hResolSubACsing=0x0;
- if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
- evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
- hResolSubBCsing=(TH1F*)lst->FindObject(hisnameBC.Data());
- cout<<hResolSubBCsing->GetName()<<endl;
- hResolSubACsing=(TH1F*)lst->FindObject(hisnameAC.Data());
- cout<<hResolSubACsing->GetName()<<endl;
- }
- Double_t error;
- Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubABsing,hResolSubBCsing,hResolSubACsing);
- Double_t resolFullmin=resolFull-error;
- Double_t resolFullmax=resolFull+error;
-
- Double_t binCentr=(Double_t)iHisC/10.+binHalfWid;
- grSingle->SetPoint(iPt,binCentr,resolFull);
- grSingle->SetPointEXlow(iPt,binHalfWid);
- grSingle->SetPointEXhigh(iPt,binHalfWid);
- grSingle->SetPointEYlow(iPt,resolFullmax-resolFull);
- grSingle->SetPointEYhigh(iPt,resolFull-resolFullmin);
- ++iPt;
- if(resolFullmin<xmin) xmin=resolFullmin;
- if(resolFullmax>xmax) xmax=resolFullmax;
- if(iHisC==minCentTimesTen){
- hResolSubAB=(TH1F*)hResolSubABsing->Clone("hResolSubAB");
- if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
- evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
- hResolSubBC=(TH1F*)hResolSubBCsing->Clone("hResolSubAB");
- hResolSubAC=(TH1F*)hResolSubACsing->Clone("hResolSubAB");
- }
- }else{
- hResolSubAB->Add(hResolSubABsing);
- if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
- evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
- hResolSubBC->Add(hResolSubBCsing);
- hResolSubAC->Add(hResolSubACsing);
- }
- }
- printf("Adding histogram %s\n",hResolSubAB->GetName());
- }
- rcflow=xmin;
- rcfhigh=xmax;
-
-
- Double_t error;
- Double_t resolFull=ComputeEventPlaneResolution(error,hResolSubAB,hResolSubBC,hResolSubAC);
-
- grInteg->SetPoint(0,40.,resolFull);
- grInteg->SetPointEXlow(0,10);
- grInteg->SetPointEXhigh(0,10.);
- grInteg->SetPointEYlow(0,error);
- grInteg->SetPointEYhigh(0,error);
-
- TCanvas* cEP=new TCanvas("cEP","EvPlaneResol");
- cEP->Divide(1,2);
- cEP->cd(1);
- hResolSubAB->Draw();
- if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
- evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap){
- hResolSubBC->SetLineColor(2);
- hResolSubBC->Draw("same");
- hResolSubAC->SetLineColor(4);
- hResolSubAC->Draw("same");
- }
- TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
- tres->SetNDC();
- tres->Draw();
- cEP->cd(2);
- grSingle->SetMarkerStyle(20);
- grInteg->SetMarkerColor(kRed+1);
- grInteg->SetLineColor(kRed+1);
- grInteg->SetMarkerStyle(29);
- grSingle->Draw("AP");
- grSingle->GetXaxis()->SetTitle("Centrality Percentile");
- grSingle->GetYaxis()->SetTitle("Resolution Correction Factor");
- grInteg->Draw("Psame");
-
- return resolFull;
-
-}
-
-//______________________________________________________________________________
-Double_t ComputeEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
- Double_t resolFull;
- if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub ||
- evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){
- resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
- error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
- }else if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
- if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta){
- resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
- error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
- }else if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
- resolFull=AliVertexingHFUtils::GetSubEvResol(hsubev1);
- error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hsubev1));
- }
- }else{
- Double_t resolSub[3];
- Double_t errors[3];
- TH1F* hevplresos[3];
- hevplresos[0]=hsubev1;
- hevplresos[1]=hsubev2;
- hevplresos[2]=hsubev3;
- for(Int_t ires=0;ires<3;ires++){
- resolSub[ires]=hevplresos[ires]->GetMean();
- errors[ires]=hevplresos[ires]->GetMeanError();
- }
- Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
- if(evPlane==AliEventPlaneResolutionHandler::kVZEROC ||
- evPlane==AliEventPlaneResolutionHandler::kVZERO){
- resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
- error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
- }
- else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){
- resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
- error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
- }
- else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta ||
- evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
- resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
- error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
- }
- }
- return resolFull;
-}
-
-//______________________________________________________________________________
-void LoadMassHistos(TList* lst, TH2F** hMassDphi){
-
- Int_t minCentTimesTen=minCent*10;
- Int_t maxCentTimesTen=maxCent*10;
- for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin]; iPtBin++){
- TString hisname=Form("hMc2deltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
- TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
- if(hMassDphi[iFinalPtBin]==0x0){
- hMassDphi[iFinalPtBin]=(TH2F*)htmp->Clone(Form("hMassCos2DphiBin%d",iFinalPtBin));
- }else{
- hMassDphi[iFinalPtBin]->Add(htmp);
- }
- printf("Adding histogram %s to PtBin %d\n",hisname.Data(),iFinalPtBin);
- }
- }
- }
-}
-
-//______________________________________________________________________________
-Bool_t DefinePtBins(TDirectoryFile* df){
-
- AliRDHFCuts *d0cuts;
- if(partname.Contains("Dzero")) d0cuts=(AliRDHFCutsD0toKpi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
- if(partname.Contains("Dplus")) d0cuts=(AliRDHFCutsDplustoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName());
- if(partname.Contains("Dstar")) d0cuts=((AliRDHFCutsDStartoKpipi*)df->Get(df->GetListOfKeys()->At(2)->GetName()));
-
- Int_t nPtBinsCuts=d0cuts->GetNPtBins();
- printf("Number of pt bins for cut object = %d\n",nPtBinsCuts);
- Float_t *ptlimsCuts=d0cuts->GetPtBinLimits();
- for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
- for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[iFinalPtBin])<0.0001){
- minPtBin[iFinalPtBin]=iPtCuts;
- if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
- }
- }
- if(TMath::Abs(ptlimsCuts[iPtCuts]-ptlims[nFinalPtBins])<0.0001) maxPtBin[nFinalPtBins-1]=iPtCuts-1;
- }
- if(TMath::Abs(ptlims[nFinalPtBins]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nFinalPtBins-1]=nPtBinsCuts-1;
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
- if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
- }
-
- return kTRUE;
-}
-
-
-//______________________________________________________________________________
-void SystForSideBands(){
- // Compute systematic ucnertainty for the side band subtraction method
- gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
-
- TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
- TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
-
- TFile *f = TFile::Open(filename.Data());
- if(!f){
- printf("file %s not found, please check file name\n",filename.Data());
- return;
- }
- TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
- if(!dir){
- printf("Directory %s not found, please check dir name\n",dirname.Data());
- return;
- }
- Bool_t binOK=DefinePtBins(dir);
- if(!binOK){
- printf("ERROR: mismatch in pt binning\n");
- return;
- }
- TList *lst =(TList*)dir->Get(listname.Data());
- if(!lst){
- printf("list %s not found in file, please check list name\n",listname.Data());
- return;
- }
-
- AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
- epres->SetEventPlane(evPlane);
- epres->SetResolutionOption(evPlaneRes);
- epres->SetUseNcollWeights(useNcollWeight);
- Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
- delete epres;
- printf("Event plane resolution %f\n",resolFull);
-
- TH2F** hMassDphi=new TH2F*[nFinalPtBins];
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- hMassDphi[iFinalPtBin]=0x0;
- }
- LoadMassHistos(lst, hMassDphi);
-
- Int_t nSteps=21;
-
- TGraphErrors** gSystSigRange=new TGraphErrors*[nFinalPtBins];
- TGraphErrors** gSystBkgRange=new TGraphErrors*[nFinalPtBins];
- TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
- TGraphErrors** gSystWhichSide=new TGraphErrors*[nFinalPtBins];
-
- Double_t min1[nFinalPtBins],max1[nFinalPtBins];
- Double_t min123[nFinalPtBins],max123[nFinalPtBins];
- Double_t min2[nFinalPtBins],max2[nFinalPtBins];
- Double_t min3[nFinalPtBins],max3[nFinalPtBins];
- Double_t min4[nFinalPtBins],max4[nFinalPtBins];
-
- Double_t sigmaRangeForSigOrig=sigmaRangeForSig;
- Double_t sigmaRangeForBkgOrig=sigmaRangeForBkg;
-
-
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("**************** BIN %d ******************\n",iFinalPtBin);
-
- Int_t rebinHistoSideBandsOrig=rebinHistoSideBands[iFinalPtBin];
-
- gSystSigRange[iFinalPtBin]=new TGraphErrors(0);
- gSystSigRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Signal Region Ptbin %d",iFinalPtBin));
-
- TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
-
- Bool_t out=FitMassSpectrum(hMass,0x0);
- if(!out) continue;
-
- min1[iFinalPtBin]=99999.;
- max1[iFinalPtBin]=-99999.;
- min123[iFinalPtBin]=99999.;
- max123[iFinalPtBin]=-99999.;
- for(Int_t iStep=0; iStep<nSteps; iStep++){
- Int_t index=iFinalPtBin*nSteps+iStep;
- sigmaRangeForSig=1.5+0.1*iStep;
- sigmaRangeForBkg=sigmaRangeForBkgOrig;
- rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
- TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
- Double_t v2obsM1=hCos2PhiSig->GetMean();
- Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
- delete hCos2PhiSig;
- Double_t v2M1=v2obsM1/resolFull;
- Double_t errv2M1=errv2obsM1/resolFull;
- gSystSigRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForSig,v2M1);
- gSystSigRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
- if(v2M1>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M1;
- if(v2M1<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M1;
- if(sigmaRangeForSig>=2. && sigmaRangeForSig<=3){
- if(v2M1>max123[iFinalPtBin]) max123[iFinalPtBin]=v2M1;
- if(v2M1<min123[iFinalPtBin]) min123[iFinalPtBin]=v2M1;
- }
- }
-
- min2[iFinalPtBin]=99999.;
- max2[iFinalPtBin]=-99999.;
- gSystBkgRange[iFinalPtBin]=new TGraphErrors(0);
- gSystBkgRange[iFinalPtBin]->SetTitle(Form("v2 vs. nSigma Bkg Region Ptbin %d",iFinalPtBin));
- for(Int_t iStep=0; iStep<nSteps; iStep++){
- Int_t index=nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
- sigmaRangeForSig=sigmaRangeForSigOrig;
- sigmaRangeForBkg=4.+0.1*iStep;
- rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
- TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
- Double_t v2obsM1=hCos2PhiSig->GetMean();
- Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
- delete hCos2PhiSig;
- Double_t v2M1=v2obsM1/resolFull;
- Double_t errv2M1=errv2obsM1/resolFull;
- gSystBkgRange[iFinalPtBin]->SetPoint(iStep,sigmaRangeForBkg,v2M1);
- gSystBkgRange[iFinalPtBin]->SetPointError(iStep,0.,errv2M1);
- if(v2M1>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M1;
- if(v2M1<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M1;
- }
-
- min3[iFinalPtBin]=99999.;
- max3[iFinalPtBin]=-99999.;
- gSystRebin[iFinalPtBin]=new TGraphErrors(0);
- gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin Ptbin %d",iFinalPtBin));
- Int_t nPts=0;
- for(Int_t iStep=0; iStep<nSteps; iStep++){
- Int_t index=2*nSteps*nFinalPtBins+iFinalPtBin*nSteps+iStep;
- sigmaRangeForSig=sigmaRangeForSigOrig;
- sigmaRangeForBkg=sigmaRangeForBkgOrig;
- rebinHistoSideBands[iFinalPtBin]=iStep+1;
- if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistoSideBands[iFinalPtBin]!=0) continue;
- TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0);
- Double_t v2obsM1=hCos2PhiSig->GetMean();
- Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
- delete hCos2PhiSig;
- Double_t v2M1=v2obsM1/resolFull;
- Double_t errv2M1=errv2obsM1/resolFull;
- gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistoSideBands[iFinalPtBin],v2M1);
- gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
- if(v2M1>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M1;
- if(v2M1<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M1;
- ++nPts;
- }
-
- min4[iFinalPtBin]=99999.;
- max4[iFinalPtBin]=-99999.;
- gSystWhichSide[iFinalPtBin]=new TGraphErrors(0);
- gSystWhichSide[iFinalPtBin]->SetTitle(Form("v2 vs. WhichSide Ptbin %d",iFinalPtBin));
- nPts=0;
- for(Int_t iStep=-1; iStep<=1; iStep++){
- Int_t index=3*nSteps*nFinalPtBins+2+iStep;
- sigmaRangeForSig=sigmaRangeForSigOrig;
- sigmaRangeForBkg=sigmaRangeForBkgOrig;
- rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
- TH1F* hCos2PhiSig=DoSideBands(index,hMassDphi[iFinalPtBin],hMass,rebinHistoSideBands[iFinalPtBin],0x0,iStep);
- Double_t v2obsM1=hCos2PhiSig->GetMean();
- Double_t errv2obsM1=hCos2PhiSig->GetMeanError();
- delete hCos2PhiSig;
- Double_t v2M1=v2obsM1/resolFull;
- Double_t errv2M1=errv2obsM1/resolFull;
- gSystWhichSide[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M1);
- gSystWhichSide[iFinalPtBin]->SetPointError(nPts,0.,errv2M1);
- if(v2M1>max4[iFinalPtBin]) max4[iFinalPtBin]=v2M1;
- if(v2M1<min4[iFinalPtBin]) min4[iFinalPtBin]=v2M1;
- ++nPts;
- }
- rebinHistoSideBands[iFinalPtBin]=rebinHistoSideBandsOrig;
- }
-
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("------ Pt Bin %d ------\n",iFinalPtBin);
- printf("Range of values for sig variation = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
- printf(" (limited to 2-3 sigma) = %f %f\n",min123[iFinalPtBin],max123[iFinalPtBin]);
- printf("Range of values for bkg variation = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
- printf("Range of values for rebin = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
- printf("Range of values for whichside = %f %f\n",min4[iFinalPtBin],max4[iFinalPtBin]);
- Float_t minenv=TMath::Min(min123[iFinalPtBin],TMath::Min(min2[iFinalPtBin],min4[iFinalPtBin]));
- Float_t maxenv=TMath::Max(max123[iFinalPtBin],TMath::Max(max2[iFinalPtBin],max4[iFinalPtBin]));
- printf(" ---> Envelope = %f %f\n",minenv,maxenv);
- }
-
- TCanvas* cs1=new TCanvas("cs1");
- cs1->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs1->cd(iFinalPtBin+1);
- gSystSigRange[iFinalPtBin]->SetMarkerStyle(20);
- gSystSigRange[iFinalPtBin]->Draw("AP");
- gSystSigRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaSignal");
- gSystSigRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-
- TCanvas* cs2=new TCanvas("cs2");
- cs2->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs2->cd(iFinalPtBin+1);
- gSystBkgRange[iFinalPtBin]->SetMarkerStyle(20);
- gSystBkgRange[iFinalPtBin]->Draw("AP");
- gSystBkgRange[iFinalPtBin]->GetXaxis()->SetTitle("nSigmaBackground");
- gSystBkgRange[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-
- TCanvas* cs3=new TCanvas("cs3");
- cs3->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs3->cd(iFinalPtBin+1);
- gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
- gSystRebin[iFinalPtBin]->Draw("AP");
- gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
- gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-
- TCanvas* cs4=new TCanvas("cs4");
- cs4->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs4->cd(iFinalPtBin+1);
- gSystWhichSide[iFinalPtBin]->SetMarkerStyle(20);
- gSystWhichSide[iFinalPtBin]->Draw("AP");
- gSystWhichSide[iFinalPtBin]->GetXaxis()->SetTitle("Side band used");
- gSystWhichSide[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-}
-
-//______________________________________________________________________________
-void SystForFitv2Mass(){
- // Compute systematic ucnertainty for the v2(M) method
- gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
-
- TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
- TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
-
- TFile *f = TFile::Open(filename.Data());
- if(!f){
- printf("file %s not found, please check file name\n",filename.Data());
- return;
- }
- TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
- if(!dir){
- printf("Directory %s not found, please check dir name\n",dirname.Data());
- return;
- }
- Bool_t binOK=DefinePtBins(dir);
- if(!binOK){
- printf("ERROR: mismatch in pt binning\n");
- return;
- }
- TList *lst =(TList*)dir->Get(listname.Data());
- if(!lst){
- printf("list %s not found in file, please check list name\n",listname.Data());
- return;
- }
-
- AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
- epres->SetEventPlane(evPlane);
- epres->SetResolutionOption(evPlaneRes);
- epres->SetUseNcollWeights(useNcollWeight);
- Double_t resolFull=epres->GetEventPlaneResolution(minCent,maxCent);
- delete epres;
- printf("Event plane resolution %f\n",resolFull);
-
- TH2F** hMassDphi=new TH2F*[nFinalPtBins];
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- hMassDphi[iFinalPtBin]=0x0;
- }
- LoadMassHistos(lst, hMassDphi);
-
- Int_t nSteps=11;
-
- TGraphErrors** gSystParamErr=new TGraphErrors*[nFinalPtBins];
- TGraphErrors** gSystRebin=new TGraphErrors*[nFinalPtBins];
- TGraphErrors** gSystLinConst=new TGraphErrors*[nFinalPtBins];
-
- Double_t min1[nFinalPtBins],max1[nFinalPtBins];
- Double_t min2[nFinalPtBins],max2[nFinalPtBins];
- Double_t min3[nFinalPtBins],max3[nFinalPtBins];
-
-
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("**************** BIN %d ******************\n",iFinalPtBin);
-
- Int_t rebinHistov2MassOrig=rebinHistov2Mass[iFinalPtBin];
- TH1F* hMass=(TH1F*)hMassDphi[iFinalPtBin]->ProjectionY();
-
- Bool_t out=FitMassSpectrum(hMass,0x0);
- if(!out) continue;
-
- min1[iFinalPtBin]=99999.;
- max1[iFinalPtBin]=-99999.;
- gSystRebin[iFinalPtBin]=new TGraphErrors(0);
- gSystRebin[iFinalPtBin]->SetTitle(Form("v2 vs. Rebin PtBin %d",iFinalPtBin));
- Int_t nPts=0;
- for(Int_t iStep=0; iStep<nSteps; iStep++){
- Int_t index=iStep;
- rebinHistov2Mass[iFinalPtBin]=iStep+1;
- if((hMassDphi[iFinalPtBin]->GetNbinsY())%rebinHistov2Mass[iFinalPtBin]!=0) continue;
- TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,useConstantvsBkgVsMass);
- Double_t v2obsM2=fv2->GetParameter(3);
- Double_t errv2obsM2=fv2->GetParError(3);
- delete fv2;
- Double_t v2M2=v2obsM2/resolFull;
- Double_t errv2M2=errv2obsM2/resolFull;
- gSystRebin[iFinalPtBin]->SetPoint(nPts,(Double_t)rebinHistov2Mass[iFinalPtBin],v2M2);
- gSystRebin[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
- if(v2M2>max1[iFinalPtBin]) max1[iFinalPtBin]=v2M2;
- if(v2M2<min1[iFinalPtBin]) min1[iFinalPtBin]=v2M2;
- ++nPts;
- }
- rebinHistov2Mass[iFinalPtBin]=rebinHistov2MassOrig;
-
- min2[iFinalPtBin]=99999.;
- max2[iFinalPtBin]=-99999.;
- gSystParamErr[iFinalPtBin]=new TGraphErrors(0);
- gSystParamErr[iFinalPtBin]->SetTitle(Form("v2 vs. ParamErr PtBin %d",iFinalPtBin));
- nPts=0;
- for(Int_t iStep=-1; iStep<=1; iStep++){
- Int_t index=nSteps*2+iStep;
- TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,iStep,useConstantvsBkgVsMass);
- Double_t v2obsM2=fv2->GetParameter(3);
- Double_t errv2obsM2=fv2->GetParError(3);
- delete fv2;
- Double_t v2M2=v2obsM2/resolFull;
- Double_t errv2M2=errv2obsM2/resolFull;
- gSystParamErr[iFinalPtBin]->SetPoint(nPts,(Double_t)iStep,v2M2);
- gSystParamErr[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
- if(v2M2>max2[iFinalPtBin]) max2[iFinalPtBin]=v2M2;
- if(v2M2<min2[iFinalPtBin]) min2[iFinalPtBin]=v2M2;
- ++nPts;
- }
-
- min3[iFinalPtBin]=99999.;
- max3[iFinalPtBin]=-99999.;
- gSystLinConst[iFinalPtBin]=new TGraphErrors(0);
- gSystLinConst[iFinalPtBin]->SetTitle(Form("v2 LinVsConst Bin%d",iFinalPtBin));
- nPts=0;
- for(Int_t iStep=0; iStep<=1; iStep++){
- Int_t index=nSteps*3+iStep;
- TF1* fv2=DoFitv2vsMass(index,hMassDphi[iFinalPtBin],hMass,rebinHistov2Mass[iFinalPtBin],0x0,0,iStep);
- Double_t v2obsM2=fv2->GetParameter(3);
- Double_t errv2obsM2=fv2->GetParError(3);
- delete fv2;
- Double_t v2M2=v2obsM2/resolFull;
- Double_t errv2M2=errv2obsM2/resolFull;
- gSystLinConst[iFinalPtBin]->SetPoint(nPts,1.-(Double_t)iStep,v2M2);
- gSystLinConst[iFinalPtBin]->SetPointError(nPts,0.,errv2M2);
- if(v2M2>max3[iFinalPtBin]) max3[iFinalPtBin]=v2M2;
- if(v2M2<min3[iFinalPtBin]) min3[iFinalPtBin]=v2M2;
- ++nPts;
- }
-
-
- }
-
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- printf("------ Pt Bin %d ------\n",iFinalPtBin);
- printf("Range of values for rebin = %f %f\n",min1[iFinalPtBin],max1[iFinalPtBin]);
- printf("Range of values for par err = %f %f\n",min2[iFinalPtBin],max2[iFinalPtBin]);
- printf("Range of values for lin const = %f %f\n",min3[iFinalPtBin],max3[iFinalPtBin]);
- Float_t minenv=TMath::Min(min2[iFinalPtBin],min3[iFinalPtBin]);
- Float_t maxenv=TMath::Max(max2[iFinalPtBin],max3[iFinalPtBin]);
- printf(" ---> Envelope = %f %f\n",minenv,maxenv);
- }
-
- TCanvas* cs1=new TCanvas("cs1");
- cs1->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs1->cd(iFinalPtBin+1);
- gSystRebin[iFinalPtBin]->SetMarkerStyle(20);
- gSystRebin[iFinalPtBin]->Draw("AP");
- gSystRebin[iFinalPtBin]->GetXaxis()->SetTitle("Rebin factor");
- gSystRebin[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-
- TCanvas* cs2=new TCanvas("cs2");
- cs2->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs2->cd(iFinalPtBin+1);
- gSystParamErr[iFinalPtBin]->SetMarkerStyle(20);
- gSystParamErr[iFinalPtBin]->Draw("AP");
- gSystParamErr[iFinalPtBin]->GetXaxis()->SetTitle("Error on Signal yield");
- gSystParamErr[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-
- TCanvas* cs3=new TCanvas("cs3");
- cs3->Divide(2,2);
- for(Int_t iFinalPtBin=0; iFinalPtBin<nFinalPtBins; iFinalPtBin++){
- cs3->cd(iFinalPtBin+1);
- gSystLinConst[iFinalPtBin]->SetMarkerStyle(20);
- gSystLinConst[iFinalPtBin]->Draw("AP");
- gSystLinConst[iFinalPtBin]->GetXaxis()->SetLimits(-0.1,1.1);
- gSystLinConst[iFinalPtBin]->GetXaxis()->SetTitle("Const/Linear v2 of background");
- gSystLinConst[iFinalPtBin]->GetYaxis()->SetTitle("v2");
- }
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff