#include "AliVertexingHFUtils.h"
#include "AliRDHFCutsDplustoKpipi.h"
#include "AliRDHFCutsDStartoKpipi.h"
+#include "AliEventPlaneResolutionHandler.h"
#include <AliVertexingHFUtils.h>
#endif
// Francesco Prino prino@to.infn.it
//global variables to be set
-TString filename="AnalysisResults_train60.root";
+TString filename="AnalysisResults_train63.root";
TString suffix="v2Dplus3050Cut4upcutPIDTPC";
TString partname="Dplus";
Int_t minCent=30;
Int_t maxCent=50;
-//evPlane flag: -1=V0C,0=V0,1=V0A,2=TPC2subevs,3=TPC3subevs
-Int_t evPlane=2;
+//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;
// pt and phi binning
const Int_t nptbinsnew=4;
// mass fit configuration
Int_t rebin[nptbinsnew]={4,4,4,4};
Int_t typeb=0; // Background: 0=expo, 1=linear, 2=pol2
+Bool_t fixAlsoMass=kFALSE;
Double_t minMassForFit=1.6;
Double_t maxMassForFit=2.2;
Float_t massRangeForCounting=0.1; // GeV
Bool_t gnopng=kFALSE; //don't save in png format (only root and eps)
Int_t minPtBin[nptbinsnew]={-1,-1,-1,-1};
Int_t maxPtBin[nptbinsnew]={-1,-1,-1,-1};
-Double_t effInOverEffOut=1.025;
+Double_t effInOverEffOut=1.03;
Double_t massD;
//methods
//method implementation
//_________________________________________________________________
Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
- Double_t resolFull;
- if(evPlane==2){
+ Double_t resolFull=1.;
+ 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];
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<=0){
+ 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==1){
+ 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==3){
+ 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]);
}
TGraphErrors* gResolVsCent=new TGraphErrors(0);
Int_t iPt=0;
Int_t nSubRes=1;
- if(evPlane!=2)nSubRes=3;//3 sub-events method for VZERO EP
+ if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+ evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
TString namereso[3]={"Reso","Reso2","Reso3"};
TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
TH2F* hevpls=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",minCentTimesTen,minCentTimesTen+25));
// TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));QUI!!!
TH1F* hevplresos[3];
Int_t ncBin=minCentTimesTen/25;
+
for(Int_t ires=0;ires<nSubRes;ires++){
hevplresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),minCentTimesTen,minCentTimesTen+25));
if(hevplresos[ires]){
fitter.DrawHere(cPhiInteg->cd(ipt+1),3,1);
}
Double_t sigma=fitter.GetSigma();
+ Double_t massFromFit=fitter.GetMean();
for(Int_t iphi=0;iphi<nphi;iphi++){
Int_t ipad=GetPadNumber(ipt,iphi);
TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
AliHFMassFitter fitter2(histtofit,hmin,hmax,1,typeb);
fitter2.SetInitialGaussianMean(massD);
fitter2.SetFixGaussianSigma(sigma);
+ if(fixAlsoMass) fitter2.SetFixGaussianMean(massFromFit);
Bool_t ok2=fitter2.MassFitter(kFALSE);
Double_t signal=0,esignal=0;
if(ok2){
TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
- if(evPlane>3||evPlane<-1){
- printf("Wrong EP flag %d \n",evPlane);
- return;
- }
-
AliRDHFCuts *cutsobj=0x0;
//Load input data from AliAnalysisTaskSEHFv2
TFile *f = TFile::Open(filename.Data());
if(!DefinePtBins(cutsobj)){
printf("cut not define pt bins\n");return;
}
+
//Load mass histograms corresponding to the required centrality, pt range and phi binning
printf("Load mass histos \n");
TList *histlist=LoadMassHistos(list,inoutanis);
Float_t massFromFit=fitter.GetMean();
Float_t sigmaFromFit=fitter.GetSigma();
TF1* funcB2=fitter.GetBackgroundRecalcFunc();
- utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2,2.5,4.5,1);
+ utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2,2.5,4.5,0.,3.,1);
}
printf("Average pt\n");
for(Int_t ipt=0;ipt<nptbinsnew;ipt++) printf("%f +- %f\n",averagePt[ipt],errorPt[ipt]);
FillSignalGraph(histlist,gSignal,gSignalfs,gSignalBC1,gSignalBC2,inoutanis);
//EP resolution
- TList* resolhist=LoadResolutionHistos(list);
- TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
- TH1F* hevplresos[3];
- TString namereso[3]={"Reso","Reso2","Reso3"};
- Int_t nSubRes=1;
- if(evPlane!=2)nSubRes=3;//3 sub-events method for VZERO EP
- for(Int_t ires=0;ires<nSubRes;ires++){
- hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
- }
- Double_t errorres;
- Double_t resol=GetEventPlaneResolution(errorres,hevplresos[0],hevplresos[1],hevplresos[2]);
+ // TList* resolhist=LoadResolutionHistos(list);
+ // TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
+ // TH1F* hevplresos[3];
+ // TString namereso[3]={"Reso","Reso2","Reso3"};
+ // Int_t nSubRes=1;
+ // if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+ // evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
+ // for(Int_t ires=0;ires<nSubRes;ires++){
+ // hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
+ // }
+ // Double_t errorres;
+ // Double_t resol=GetEventPlaneResolution(errorres,hevplresos[0],hevplresos[1],hevplresos[2]);
+ AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
+ epres->SetEventPlane(evPlane);
+ epres->SetResolutionOption(evPlaneRes);
+ epres->SetUseNcollWeights(useNcollWeight);
+ Double_t resol=epres->GetEventPlaneResolution(minCent,maxCent);
+ delete epres;
printf("Event plane resolution %f\n",resol);
printf("Compute v2 \n");
TList *list =(TList*)dir->Get(listname.Data());
if(!list){
printf("list %s not found in file, please check list name\n",listname.Data());return;
+ }
+ if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta &&
+ evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
+ printf("Wrong setting of event plane resolution forced it to kTwoSub\n");
+ evPlaneRes=AliEventPlaneResolutionHandler::kTwoRandSub;
}
TList* resolhist=LoadResolutionHistos(list);
TGraphErrors* gResolVsCent=(TGraphErrors*)resolhist->FindObject("gResolVsCent");
TH1F* hevplresos[3];
TString namereso[3]={"Reso","Reso2","Reso3"};
Int_t nSubRes=1;
- if(evPlane!=2)nSubRes=3;//3 sub-events method for VZERO EP
+ if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
+ evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
for(Int_t ires=0;ires<nSubRes;ires++){
hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
}
// gStyle->SetTitleAlign(22);
gStyle->SetTitleFont(42,"xyzg");
+ TH1F* htpc = (TH1F*)hevpls->ProjectionX();
+ TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
+ TH1F* hUsedPlane;
+ if(evPlane==AliEventPlaneResolutionHandler::kVZERO ||
+ evPlane==AliEventPlaneResolutionHandler::kVZEROA ||
+ evPlane==AliEventPlaneResolutionHandler::kVZEROC) hUsedPlane=(TH1F*)hv0->Clone("hUsedPlane");
+ else hUsedPlane=(TH1F*)htpc->Clone("hUsedPlane");
+
TCanvas* cvtotevpl=new TCanvas(Form("cvtotevpl%s",suffixcentr.Data()),Form("Ev plane %s",suffixcentr.Data()),1200,400);
cvtotevpl->Divide(3,1);
cvtotevpl->cd(1);
hevpls->SetTitleFont(42);
hevpls->Draw("COLZ");
- TH1F* htpc = (TH1F*)hevpls->ProjectionX();
cvtotevpl->cd(2);
htpc->Draw();
htpc->Fit("pol0");
- TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
cvtotevpl->cd(3);
hv0->Draw();
hv0->Fit("pol0");
+ TCanvas* cvfitevpl=new TCanvas(Form("cvditevpl%s",suffixcentr.Data()),Form("Fit to Ev plane %s",suffixcentr.Data()));
+ cvfitevpl->cd();
+ hUsedPlane->Draw();
+ TF1* four2=new TF1("four2","[0]*(1+2*[1]*cos(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
+ TF1* four2s=new TF1("four2s","[0]*(1+2*[1]*sin(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
+ hUsedPlane->Fit(four2s);
+ hUsedPlane->Fit(four2);
+ four2->SetLineColor(1);
+ four2s->SetLineColor(4);
+ four2->Draw("same");
+ four2s->Draw("same");
+ hUsedPlane->SetMaximum(hUsedPlane->GetMaximum()*1.07);
+ TLatex* tsin=new TLatex(0.15,0.84,Form("1+2*(%.4f)*sin(2*#Phi_{EP})",four2s->GetParameter(1)));
+ tsin->SetNDC();
+ tsin->SetTextColor(4);
+ tsin->Draw();
+ TLatex* tcos=new TLatex(0.15,0.77,Form("1+2*(%.4f)*cos(2*#Phi_{EP})",four2->GetParameter(1)));
+ tcos->SetNDC();
+ tcos->SetTextColor(1);
+ tcos->Draw();
+
+ // Compute the second Fourier component for since and cosine
+ Double_t aveCos2Phi=0.;
+ Double_t aveSin2Phi=0.;
+ Double_t counts=0.;
+ for(Int_t i=1; i<=hUsedPlane->GetNbinsX(); i++){
+ Double_t co=hUsedPlane->GetBinContent(i);
+ Double_t phi=hUsedPlane->GetBinCenter(i);
+ aveCos2Phi+=co*TMath::Cos(2*phi);
+ aveSin2Phi+=co*TMath::Sin(2*phi);
+ counts+=co;
+ }
+ if(counts>0){
+ aveCos2Phi/=counts;
+ aveSin2Phi/=counts;
+ }
+ printf("\n------ Fourier 2nd components of EP distribution ------\n");
+ printf("<cos(2*Psi_EP)>=%f\n",aveCos2Phi);
+ printf("<sin(2*Psi_EP)>=%f\n",aveSin2Phi);
+ printf("\n");
+
+
TCanvas* cvtotevplreso=new TCanvas(Form("cvtotevplreso%s",suffixcentr.Data()),Form("Ev plane Resolution %s",suffixcentr.Data()));
cvtotevplreso->cd();
hevplresos[0]->SetTitle("TPC cos2(#Psi_{A}-#Psi_{B})");
Double_t error;
Double_t resolFull=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
- TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+ TPaveText* pvreso=new TPaveText(0.1,0.35,0.6,0.48,"NDC");
pvreso->SetBorderSize(0);
pvreso->SetFillStyle(0);
+ pvreso->AddText(Form("Number of events = %.0f\n",hevplresos[0]->GetEntries()));
pvreso->AddText(Form("Resolution on full event = %.4f#pm%.4f\n",resolFull,error));
pvreso->Draw();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff