PWG3/vertexingHF/charmFlow/DmesonsFlowAnalysis.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TFile.h>
   4 #include <TString.h>
   5 #include <TH2F.h>
   6 #include <TH1F.h>
   7 #include <TF1.h>
   8 #include <TGraph.h>
   9 #include <TGraphErrors.h>
  10 #include <TMultiGraph.h>
  11 #include <TDirectoryFile.h>
  12 #include <TList.h>
  13 #include <TCanvas.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TPaveText.h>
  17 #include <TStyle.h>
  18 #include <TASImage.h>
  19 #include <TPad.h>
  20 #include <TROOT.h>
  21 #include <TDatabasePDG.h>
  22 #include <TParameter.h>
  23 #include <TLatex.h>
  24
  25 #include <AliHFMassFitter.h>
  26 #include "AliRDHFCutsD0toKpi.h"
  27 #include "AliRDHFCutsDplustoKpipi.h"
  28 #include "AliRDHFCutsDStartoKpipi.h"
  29
  30 #endif
  31
  32 //methods for the analysis of AliAnalysisTaskSEHFv2 output
  33 //Authors: Chiara Bianchin cbianchi@pd.infn.it
  34 //         Giacomo Ortona  ortona@to.infn.it
  35 //         Francesco Prino prino@to.infn.it
  36
  37 //global variables to be set
  38 // const Int_t nptbinsnew=6;
  39 // Float_t ptbinsnew[nptbinsnew+1]={2,3,4,5,8,12,999.};
  40 const Int_t nptbinsnew=3;
  41 Float_t ptbinsnew[nptbinsnew+1]={2,5,8,16.};
  42 //const Int_t nptbinsnew=4;
  43 //Float_t ptbinsnew[nptbinsnew+1]={2,4,5,8,999.};
  44 //Float_t ptbinsnew[nptbinsnew+1]={2,3,5,8,999.};
  45 // const Int_t nptbinsnew=5;
  46 // Float_t ptbinsnew[nptbinsnew+1]={2,3,4,5,8,999.};
  47 // const Int_t nptbinsnew=2;
  48 // Float_t ptbinsnew[nptbinsnew+1]={2,8,999.};
  49 Int_t fittype=0;
  50 //Int_t rebin[nptbinsnew]={4,4,4,5};
  51 Int_t rebin[nptbinsnew]={4,4,4};
  52 Double_t nsigma=3;
  53
  54 //methods
  55 Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutobj,TString listname,TString partname,TString path="./",TString filename="AnalysisResults.root");
  56 void WriteCanvas(TCanvas* cv,TString text);
  57 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
  58 void InOutPic(TVirtualPad *c,Int_t inout=0,TString where="tr");//inout: 0=IN, 1=OUT
  59 void PhiBinPic(TVirtualPad *c,Int_t angle,TString where);
  60 Double_t FindChi(Double_t res, Int_t k);
  61 Double_t Pol(Double_t x, Int_t k);
  62 Double_t ResolK1(Double_t x);
  63 Double_t ComputeResol(Double_t resSub, Int_t k);
  64
  65
  66 //methods implementation
  67
  68 Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutsobj,TString listname,TString partname,TString path,TString filename){
  69
  70   TString hstatname="hEventsInfo",dirname="PWG3_D2H_HFv2";
  71   filename.Prepend(path);
  72   listname+=partname;
  73   hstatname+=partname;
  74   // TString tmpsuff="NoCos"; //"Nod0d0"
  75   // listname+=tmpsuff;
  76   // hstatname+=tmpsuff;
  77
  78   TFile* f=new TFile(filename.Data());
  79   if(!f){
  80     cout<<filename.Data()<<" not found"<<endl;
  81     return kFALSE;
  82   }
  83   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
  84   if(!f){
  85     cout<<dirname.Data()<<" not found  in "<<filename.Data()<<endl;
  86     return kFALSE;
  87   }
  88
  89   list=(TList*)dir->Get(listname);
  90   if(!list){
  91     cout<<"List "<<listname.Data()<<" not found"<<endl;
  92     dir->ls();
  93     return kFALSE;
  94   }
  95
  96   hstat=(TH1F*)dir->Get(hstatname);
  97   if(!hstat){
  98     cout<<hstatname.Data()<<" not found"<<endl;
  99     return kFALSE;
 100   }
 101
 102   if(partname.Contains("D0")) cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 103   if(partname.Contains("Dplus")) cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 104   if(partname.Contains("Dstar")) cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 105
 106   if(!cutsobj){
 107     cout<<"Cut object not found! check position of the key!"<<endl;
 108     return kFALSE;
 109   }
 110
 111   return kTRUE;
 112 }
 113
 114 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
 115   for (Int_t i=0;i<nbins;i++){
 116     //cout<<value<<" from "<<array[i]<<" to "<<array[i+1]<<"?"<<endl;
 117     if(value>=array[i] && value<array[i+1]){
 118       return i;
 119     }
 120   }
 121   cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
 122   return -1;
 123
 124 }
 125
 126 void WriteCanvas(TCanvas* cv,TString text){
 127   cv->SaveAs(Form("%s%s.png",cv->GetName(),text.Data()));
 128 }
 129
 130 void InOutPic(TVirtualPad *c,Int_t inout,TString where){
 131   TPad *myPadLogo = 0x0;
 132   if(where=="tr"){
 133     myPadLogo = new TPad("myPadPic", "Pad for In/Out pic",0.66,0.68,0.86,0.88);
 134   }
 135
 136   myPadLogo->SetFillColor(0);
 137   myPadLogo->SetBorderMode(0);
 138   myPadLogo->SetBorderSize(2);
 139   myPadLogo->SetFrameBorderMode(0);
 140   myPadLogo->SetLeftMargin(0.0);
 141   myPadLogo->SetTopMargin(0.0);
 142   myPadLogo->SetBottomMargin(0.0);
 143   myPadLogo->SetRightMargin(0.0);
 144   myPadLogo->Draw();
 145   myPadLogo->cd();
 146   TASImage *myAliceLogo = new TASImage(Form("/home/cbianchi/macros/sketch%s.png",inout==0 ? "IN" : "OUT"));
 147   myAliceLogo->Draw();
 148   c->cd();
 149   TPaveText* t1=0x0;//=new TPaveText(0.62,0.63,0.89,0.71,"NDC");
 150   if(where=="tr"){
 151     t1=new TPaveText(0.65,0.58,0.85,0.67,"NDC");
 152   }
 153   t1->SetFillStyle(0);
 154   t1->AddText(Form("%s-PLANE",inout==0 ? "IN" : "OUT-OF"));
 155   if(inout==0)t1->SetTextColor(kBlue);
 156   else t1->SetTextColor(kRed);
 157   t1->Draw();
 158 }
 159
 160 void PhiBinPic(TVirtualPad *c,Int_t angle,TString where){
 161
 162   TString picname="angles";
 163   if(angle==0) picname.Append("0pi4");
 164   if(angle==1) picname.Append("pi4pi2");
 165   if(angle==2) picname.Append("pi232pi");
 166   if(angle==3) picname.Append("32pipi");
 167   picname.Append(".png");
 168
 169   TPad *myPadLogo = 0x0;
 170   if(where=="tr"){
 171     myPadLogo = new TPad("myPadPic", "Pad for bin of phi pic",0.5,0.75,0.8,0.8);
 172   }
 173
 174   myPadLogo->SetFillColor(0);
 175   myPadLogo->SetBorderMode(0);
 176   myPadLogo->SetBorderSize(2);
 177   myPadLogo->SetFrameBorderMode(0);
 178   myPadLogo->SetLeftMargin(0.0);
 179   myPadLogo->SetTopMargin(0.0);
 180   myPadLogo->SetBottomMargin(0.0);
 181   myPadLogo->SetRightMargin(0.0);
 182   myPadLogo->Draw();
 183   myPadLogo->cd();
 184   TASImage *myAliceLogo = new TASImage(Form("/home/cbianchi/macros/%s",picname.Data()));
 185   c->cd();
 186   myAliceLogo->Draw();
 187
 188
 189   /*
 190   TPaveText* t1=0x0;//=new TPaveText(0.62,0.63,0.89,0.71,"NDC");
 191   if(where=="tr"){
 192     t1=new TPaveText(0.65,0.58,0.85,0.67,"NDC");
 193   }
 194   t1->SetFillStyle(0);
 195   t1->AddText(Form("%s-PLANE",inout==0 ? "IN" : "OUT-OF"));
 196   if(inout==0)t1->SetTextColor(kBlue);
 197   else t1->SetTextColor(kRed);
 198   t1->Draw();
 199   */
 200 }
 201
 202
 203 void DmesonsSignalExtraction(TString partname="D0",Int_t mincentr=30,Int_t maxcentr=70,TString textleg="",TString path="./",Bool_t official=kFALSE){
 204
 205   //read the output of HFv2task and extract the signal in pt bins in plane and out of plane and in bins of phi pt integrated
 206
 207   gStyle->SetCanvasColor(0);
 208   gStyle->SetTitleFillColor(0);
 209   gStyle->SetStatColor(0);
 210   gStyle->SetPalette(1);
 211   gStyle->SetOptStat(0);
 212
 213   Int_t info=1;
 214   if (official) info=0;
 215
 216   if(textleg==""){
 217     textleg=Form("centr%d-%d",mincentr,maxcentr);
 218   }
 219   //  Double_t pi=TMath::Pi();
 220   TString listname="coutputv2";
 221
 222   TList* list;
 223   TH1F * hstat;
 224   AliRDHFCuts* cutobj;
 225
 226   Bool_t isRead=ReadFile(list,hstat,cutobj,listname,partname,path);
 227   if(!isRead) return;
 228   if(!list || !hstat || !cutobj){
 229     cout<<":-( null pointers..."<<endl;
 230     return;
 231   }
 232   //uncomment when problems with the cut obj
 233   /*
 234   const Int_t nptbins=12;
 235   Float_t ptbins[nptbins+1]={0,1,2,3,4,5,6,8,12,16,999};
 236   */
 237
 238   Float_t* ptbins=cutobj->GetPtBinLimits();
 239   const Int_t nptbins=cutobj->GetNPtBins();
 240
 241   TH1F* hnphibins=(TH1F*)list->FindObject("hPhiBins");
 242   Int_t nphibins=hnphibins->GetNbinsX();
 243   Double_t phibins[nphibins],dx[nphibins];
 244   cout<<"(";
 245   for(Int_t i=0;i<nphibins;i++){
 246     Double_t width=hnphibins->GetBinWidth(i+1);
 247     phibins[i]=hnphibins->GetBinLowEdge(i+1)+width/2.;
 248     dx[i]=width/2.;
 249     cout<<phibins[i]<<",";
 250   }
 251   cout<<")"<<endl;
 252   cout<<"Number of phi bins = "<<nphibins<<endl;
 253
 254   cout<<"Number of pt bins = "<<cutobj->GetNPtBins()<<endl;
 255   Double_t deltapt[nptbinsnew],pt[nptbinsnew],ptledges[nptbinsnew+1];
 256   for(Int_t j=0;j<nptbinsnew;j++){
 257     deltapt[j]=(ptbinsnew[j+1]-ptbinsnew[j])/2.;
 258     pt[j]=ptbinsnew[j]+deltapt[j];
 259     ptledges[j]=ptbinsnew[j];
 260
 261     if(j==nptbinsnew-1 && deltapt[j]>10) {
 262       pt[j]=ptbinsnew[j]+10;
 263       deltapt[j]=10;
 264     }
 265   }
 266   ptledges[nptbinsnew]=ptledges[nptbinsnew-1]+deltapt[nptbinsnew-1];
 267
 268   //invariant mass pt integrated, in bins of phi
 269   TH1F** hmassptint=new TH1F*[nphibins];
 270   //invariant mass in bins of pt and in/out plane
 271   TH1F*** hmassnonintphi=new TH1F**[2];
 272
 273   //init
 274   for(Int_t j=0;j<2;j++){
 275     hmassnonintphi[j]=new TH1F*[nptbinsnew];
 276     for(Int_t i=0;i<nptbinsnew;i++){
 277       hmassnonintphi[j][i]=0x0;
 278       cout<<"Init "<<hmassnonintphi[j][i]<<endl;
 279     }
 280   }
 281   for(Int_t j=0;j<nphibins;j++){
 282     hmassptint[j]=0x0;
 283   }
 284
 285
 286   //output file
 287   TFile* fout=new TFile(Form("HistoInputv2Calc%s.root",textleg.Data()),"recreate");
 288   TClonesArray ptblim("TParameter<float>",nptbinsnew+2); //the first is the number of bins
 289   new(ptblim[0])TParameter<float>("nptbins",nptbinsnew);
 290   for(Int_t ival=0;ival<nptbinsnew+1;ival++){
 291     TString name=Form("ptbin%d",ival);
 292     new(ptblim[ival+1])TParameter<float>(name.Data(),ptbinsnew[ival]);
 293   }
 294   fout->cd();
 295   ptblim.Write();
 296   cutobj->Write();
 297
 298
 299   //reading the hostograms and filling hmassnonintphi and hmassptint
 300
 301   Bool_t isInPlane=kFALSE, isOutOfPlane=kFALSE;
 302   Int_t indx=-1;
 303
 304   for(Int_t i=0;i<nphibins;i++){ //loop on phi
 305     cout<<" -- phi bin "<<i<<"/"<<nphibins-1<<endl;
 306
 307     for(Int_t j=0;j<nptbins;j++){ //loop on pt
 308       cout<<"   ** pt bin "<<j<<"/"<<nptbins-1<<endl;
 309
 310       isInPlane=kFALSE, isOutOfPlane=kFALSE;
 311
 312       TH1F* h=(TH1F*)list->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",j,i,mincentr,mincentr+5));
 313       if(!h){
 314         cout<<"qui Histogram pt "<<j<<", phi "<<i<<", centr "<<mincentr<<"-"<<mincentr+5<<" not found"<<endl;
 315         //list->ls();
 316         continue;
 317         //return;
 318       }
 319       for(Int_t icentr=mincentr+5;icentr<maxcentr;icentr=icentr+5){
 320         TH1F* hcentr05=(TH1F*)list->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",j,i,icentr,icentr+5));
 321         if(!hcentr05){
 322           cout<<"Histogram pt "<<j<<", phi "<<i<<", centr "<<icentr<<"-"<<icentr+5<<" not found"<<endl;
 323           return;
 324         }
 325         h->Add(hcentr05);
 326       }
 327       Int_t kptnew=FindPtBin(nptbinsnew,ptbinsnew,ptbins[j]);
 328       cout<<"----*** pt"<<kptnew<<endl;
 329       if(kptnew==-1) continue;
 330       //Filling in plane and out of plane
 331       if((phibins[i]>0 && phibins[i]<= 0.79) || (phibins[i]> 2.36 && phibins[i]<3.14)) {isInPlane=kTRUE; indx=0;}
 332       if(phibins[i]> 0.79 && phibins[i]<= 2.36) {isOutOfPlane=kTRUE; indx=1;cout<<" IS OUT of plane"<<endl;}
 333       if(!hmassnonintphi[indx][kptnew]) {
 334
 335         cout<<"\tpt"<<j<<", phi"<<i<<" in phi "<<i<<" ptnew "<<kptnew<<endl;
 336         hmassnonintphi[indx][kptnew]=(TH1F*)h->Clone(Form("hMass_ptnw%dphi%d",kptnew,isInPlane? 0 : 1));
 337         hmassnonintphi[indx][kptnew]->SetTitle(Form("Mass ptnw%d phi%s",kptnew,isInPlane ? "In plane" : "Out of plane"));
 338
 339       }
 340       else {
 341         cout<<"\tadding pt"<<j<<", phi"<<i<<" in phi "<<i<<" ptnew "<<kptnew<<endl;
 342         hmassnonintphi[indx][kptnew]->Add(h);
 343       }
 344
 345
 346       //pt integrated mass plots
 347       if(!hmassptint[i]) {
 348         cout<<"Filling histo ptint["<<i<<"] with hMass_pt"<<j<<"phi"<<i<<endl;
 349         hmassptint[i]=(TH1F*)h->Clone(Form("hphi_%d",i));
 350         hmassptint[i]->SetTitle(Form("Invariant Mass (Phi %d, Pt integr)",i));
 351         //cout<<"Entries "<<hmassptint[i]->GetEntries()<<endl;
 352       }
 353       else {
 354         cout<<"Adding histo hMass_pt"<<j<<"phi"<<i<<" to ptint["<<i<<"]"<<endl;
 355         hmassptint[i]->Add(h);
 356       }
 357     } //end pt loop
 358   } //end phi loop
 359
 360   //In-plane Out-of-plane anisotropy
 361
 362   //canvas
 363   //mass in bins of pt
 364   TCanvas* cvnewptbins=new TCanvas("cvnewptbins","Some bins in pt, in-plane/out-of-plane",1600,800);
 365   cvnewptbins->Divide(nptbinsnew,2);
 366   //chi square per ptbin
 367   TCanvas* cvchinewptbins=new TCanvas("cvchinewptbins", "Chi vs pt");
 368
 369   //legend IN-PLANE OUT-OF-PLANE
 370   TLegend* legsig=new TLegend(0.7,0.6,0.9,0.8,"");
 371   legsig->SetBorderSize(0);
 372   legsig->SetFillStyle(0);
 373
 374   TPaveText* txtoff=0x0;
 375   if(official){
 376     gStyle->SetOptTitle(0);
 377     gStyle->SetFrameBorderMode(0);
 378     txtoff=new TPaveText(0.1,0.53,0.55,0.83,"NDC");
 379     txtoff->SetBorderSize(0);
 380     txtoff->SetFillStyle(0);
 381     txtoff->AddText(Form("D^{0}#rightarrow K^{-}#pi^{+} %.0f#times10^{6} events",hstat->Integral(0,1)/1e6));
 382     txtoff->AddText("Pb-Pb @ #sqrt{s_{NN}} = 2.76TeV");
 383     txtoff->AddText(Form("Centrality %d-%d%s",mincentr,maxcentr,"%"));
 384     txtoff->SetTextColor(kCyan+3);
 385   }
 386
 387   Double_t signal[2][nptbinsnew],background[2][nptbinsnew],significance[2][nptbinsnew];
 388   Double_t errsgn[2][nptbinsnew],errbkg[2][nptbinsnew],errsignificance[2][nptbinsnew];
 389   Int_t k=0;
 390   //TString textinoutpl[2]={"0<#Delta#phi<#pi/4 && 3/4#pi<#Delta#phi<#pi","#pi/4<#Delta#phi<3/4#pi"};
 391
 392
 393   //histograms
 394   //chi square
 395   TH1F** hchivsptnew=new TH1F*[2];
 396
 397
 398   for(Int_t i=0;i<2;i++){ //in plane/out of plane
 399     printf("Drawing %s\n",i==0 ? "in plane" : "out of plane");
 400     hchivsptnew[i]=new TH1F(Form("hchivsptnew%d",i), "Chi square from fit;p_{t} (GeV/c);#chi^{2}",nptbinsnew,ptledges);
 401     hchivsptnew[i]->SetMarkerColor(i+1);
 402     hchivsptnew[i]->SetMarkerStyle(20);
 403
 404      for(Int_t j=0;j<nptbinsnew;j++){
 405       TPaveText* pvbinphipt=new TPaveText(0.1,0.8,0.7,0.9,"NDC");
 406       pvbinphipt->SetBorderSize(0);
 407       pvbinphipt->SetFillStyle(0);
 408       if(j!=nptbinsnew-1) pvbinphipt->AddText(Form("%.0f<p_{t}<%.0f GeV/c",ptbinsnew[j],ptbinsnew[j+1]));
 409       else pvbinphipt->AddText(Form("p_{t}>%.0f GeV/c",ptbinsnew[j]));
 410       k++;
 411
 412       if(!hmassnonintphi[i][j]) {
 413         printf("hist %s pt %d not found\n",i==0 ? "In-plane" : "out-of-plane",j);
 414         continue;
 415       }
 416
 417       fout->cd();
 418       hmassnonintphi[i][j]->Write();
 419
 420       cvnewptbins->cd(k);
 421
 422       //Fit mass histograms in bins of pt in-plane or out-of-plane
 423       AliHFMassFitter fitter(hmassnonintphi[i][j],hmassnonintphi[i][j]->GetBinLowEdge(1),hmassnonintphi[i][j]->GetBinLowEdge(hmassnonintphi[i][j]->GetNbinsX()+1),rebin[j],fittype);
 424
 425       Bool_t ok=fitter.MassFitter(kFALSE);
 426       if(ok) {
 427         fitter.DrawHere(cvnewptbins->cd(k),nsigma,info);
 428
 429         fitter.Signal(nsigma,signal[i][j],errsgn[i][j]);
 430         fitter.Background(nsigma,background[i][j],errbkg[i][j]);
 431         fitter.Significance(nsigma,significance[i][j],errsignificance[i][j]);
 432
 433         //fill chi square
 434         hchivsptnew[i]->SetBinContent(j,fitter.GetChiSquare());
 435
 436       }
 437       else { //fit failed
 438         hmassnonintphi[i][j]->Draw();
 439         hchivsptnew[i]->SetBinContent(j,-1);
 440       }
 441
 442       //draw mass
 443       pvbinphipt->Draw();
 444      }
 445
 446      //draw chi square
 447      cvchinewptbins->cd();
 448      if(i==0){
 449        hchivsptnew[i]->Draw("ptext");
 450        legsig->AddEntry(hchivsptnew[i],"In-plane","p");
 451      }
 452      else {
 453        hchivsptnew[i]->Draw("psamestext");
 454        legsig->AddEntry(hchivsptnew[i],"Out-on-plane","p");
 455      }
 456   }
 457
 458
 459   //write png of mass as a function of pt in plane and out of plane and chi square
 460   cvchinewptbins->cd();
 461   legsig->Draw();
 462   WriteCanvas(cvchinewptbins,textleg);
 463   WriteCanvas(cvnewptbins,Form("inout%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
 464
 465   //write canvas chi square
 466   fout->cd();
 467   cvchinewptbins->Write();
 468
 469   //graphs of yields
 470   TGraphErrors** gr=new TGraphErrors*[2];
 471   TGraphErrors** grb=new TGraphErrors*[2];
 472   TGraphErrors** grsgf=new TGraphErrors*[2];
 473   TGraph** grerr=new TGraph*[2];
 474   TGraph** grberr=new TGraph*[2];
 475   //canvas for these graphs
 476   TCanvas* cvsig=new TCanvas("cvsig","Signal vs pt");
 477   TCanvas* cvsgnf=new TCanvas("cvsgnf","Significance vs pt");
 478   TCanvas* cvbkg=new TCanvas("cvbkg","Background vs pt");
 479   TCanvas* cvsiger=new TCanvas("cvsiger","Error on Signal vs pt");
 480   TCanvas* cvbkger=new TCanvas("cvbkger","Error on Background vs pt");
 481
 482   for(Int_t i=0;i<2;i++) { //in/out
 483     gr[i]=new TGraphErrors(0);
 484     grb[i]=new TGraphErrors(0);
 485     grsgf[i]=new TGraphErrors(0);
 486     grerr[i]=new TGraph(0);
 487     grberr[i]=new TGraph(0);
 488
 489     for(Int_t j=0;j<nptbinsnew;j++){ //loop on ptbins
 490
 491       //signal
 492       gr[i]->SetPoint(j,pt[j],signal[i][j]);
 493       gr[i]->SetPointError(j,deltapt[j],errsgn[i][j]);
 494       cout<<j<<" pt "<<pt[j]<<" pm "<<deltapt[j]<<", sgn "<<signal[i][j]<<" pm "<<errsgn[i][j]<<endl;
 495
 496       //background
 497       grb[i]->SetPoint(j,pt[j],background[i][j]);
 498       grb[i]->SetPointError(j,deltapt[j],errbkg[i][j]);
 499       grsgf[i]->SetPoint(j,pt[j],significance[i][j]);
 500       grsgf[i]->SetPointError(j,deltapt[j],errsignificance[i][j]);
 501
 502       //error on signal
 503       grerr[i]->SetPoint(j,pt[j],errsgn[i][j]/signal[i][j]);
 504       grberr[i]->SetPoint(j,pt[j],errbkg[i][j]/background[i][j]);
 505     }
 506
 507     //titles etc.
 508     gr[i]->SetTitle("Signal vs p_{t}");
 509     gr[i]->SetName(Form("sig%s", i==0 ? "inplane" : "outofplane"));
 510     gr[i]->GetXaxis()->SetTitle("p_{t}");
 511     gr[i]->GetYaxis()->SetTitle("Signal");
 512     gr[i]->SetMinimum(0);
 513     gr[i]->SetMarkerStyle(20);
 514     gr[i]->SetMarkerColor(i+1);
 515
 516     grb[i]->SetTitle("Background vs p_{t}");
 517     grb[i]->SetName(Form("bkg%s", i==0 ? "inplane" : "outofplane"));
 518     grb[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 519     grb[i]->GetYaxis()->SetTitle("Background");
 520     grb[i]->SetMinimum(0);
 521     grb[i]->SetMarkerStyle(20);
 522     grb[i]->SetMarkerColor(i+1);
 523
 524     grsgf[i]->SetTitle("Significance vs p_{t};p_{t} (GeV/c);Significance");
 525     grsgf[i]->SetName(Form("sgf%s", i==0 ? "inplane" : "outofplane"));
 526     grsgf[i]->SetMinimum(0);
 527     grsgf[i]->SetMarkerStyle(20);
 528     grsgf[i]->SetMarkerColor(i+1);
 529
 530     grerr[i]->SetTitle("#epsilon_{r} Signal");
 531     grerr[i]->SetName(Form("sigerr%s", i==0 ? "inplane" : "outofplane"));
 532     grerr[i]->GetYaxis()->SetTitle("#epsilon_{r} Signal");
 533     grerr[i]->GetXaxis()->SetTitle("p_{t}");
 534     grerr[i]->SetMinimum(0);
 535     grerr[i]->SetMarkerStyle(20);
 536     grerr[i]->SetMarkerColor(i+1);
 537
 538     grberr[i]->SetTitle("#epsilon_{r} Background");
 539     grberr[i]->SetName(Form("bkgerr%s", i==0 ? "inplane" : "outofplane"));
 540     grberr[i]->GetYaxis()->SetTitle("#epsilon_{r} Background");
 541     grberr[i]->GetXaxis()->SetTitle("p_{t}");
 542     grberr[i]->SetMinimum(0);
 543     grberr[i]->SetMarkerStyle(20);
 544     grberr[i]->SetMarkerColor(i+1);
 545
 546     //draw signal, bkg, significance, errors
 547    if(i==0) {
 548       cvsig->cd();
 549       gr[i]->Draw("AP");
 550
 551       cvbkg->cd();
 552       grb[i]->Draw("AP");
 553
 554       cvsgnf->cd();
 555       grsgf[i]->Draw("AP");
 556
 557       cvsiger->cd();
 558       grerr[i]->Draw("AP");
 559
 560       cvbkger->cd();
 561       grberr[i]->Draw("AP");
 562     }
 563     else {
 564       cvsig->cd();
 565       gr[i]->Draw("P");
 566
 567       cvbkg->cd();
 568       grb[i]->Draw("P");
 569
 570       cvsgnf->cd();
 571       grsgf[i]->Draw("P");
 572
 573       cvsiger->cd();
 574       grerr[i]->Draw("P");
 575
 576       cvbkger->cd();
 577       grberr[i]->Draw("P");
 578    }
 579
 580    //write in output file
 581    fout->cd();
 582    gr[i]->Write();
 583    grb[i]->Write();
 584    grerr[i]->Write();
 585    grsgf[i]->Write();
 586   }
 587
 588   cvsig->cd();
 589   legsig->Draw();
 590   cvbkg->cd();
 591   legsig->Draw();
 592   cvsgnf->cd();
 593   legsig->Draw();
 594   cvsiger->cd();
 595   legsig->Draw();
 596   cvbkger->cd();
 597   legsig->Draw();
 598
 599   //write png
 600   WriteCanvas(cvsig,textleg);
 601   WriteCanvas(cvbkg,textleg);
 602   WriteCanvas(cvsiger,textleg);
 603   WriteCanvas(cvbkger,textleg);
 604   WriteCanvas(cvsgnf,textleg);
 605
 606   //fit again with sigma obtained from total inv mass (in-plane+out-of-plane)
 607   Double_t sigmafullmass[nptbinsnew];
 608
 609   //canvas
 610   TCanvas* cvmasstotnewptbins=new TCanvas("cvmasstotnewptbins", "Mass in some pt bins in+out",1600,400);
 611   cvmasstotnewptbins->Divide(nptbinsnew,1);
 612   TCanvas* cvnewptbins2=new TCanvas("cvnewptbins2", "Mass in some pt bins in/out of plane fixed sigma",1600,800);
 613   cvnewptbins2->Divide(nptbinsnew,2);
 614
 615   for(Int_t i=0;i<nptbinsnew;i++){//pt bins
 616
 617     //pave text
 618     TPaveText* pvbinpt=new TPaveText(0.1,0.8,0.7,0.9,"NDC");
 619     pvbinpt->SetBorderSize(0);
 620     pvbinpt->SetFillStyle(0);
 621     if(i!= nptbinsnew-1) pvbinpt->AddText(Form("%.0f<p_{t}<%.0f GeV/c",ptbinsnew[i],ptbinsnew[i+1]));
 622     else if(ptbinsnew[i+1]-ptbinsnew[i] > 10) pvbinpt->AddText(Form("p_{t}>%.0f GeV/c",ptbinsnew[i]));
 623
 624     //mass histogram integrated in phi
 625     TH1F* hmassfull=(TH1F*)hmassnonintphi[0][i]->Clone();
 626     hmassfull->Add(hmassnonintphi[1][i]);
 627
 628     AliHFMassFitter fitter(hmassfull,hmassfull->GetBinLowEdge(1),hmassfull->GetBinLowEdge(hmassfull->GetNbinsX()+1),rebin[i],fittype);
 629     Bool_t ok=fitter.MassFitter(kFALSE);
 630     if(ok){
 631       fitter.DrawHere(cvmasstotnewptbins->cd(i+1),nsigma,info);
 632       sigmafullmass[i]=fitter.GetSigma();
 633     }else{
 634       cvmasstotnewptbins->cd(i+1);
 635       hmassfull->Draw();
 636       sigmafullmass[i]=-2;
 637     }
 638     cvmasstotnewptbins->cd(i+1);
 639     pvbinpt->Draw();
 640
 641
 642     //in-plane
 643     AliHFMassFitter fitter0(hmassnonintphi[0][i],hmassnonintphi[0][i]->GetBinLowEdge(1),hmassnonintphi[0][i]->GetBinLowEdge(hmassnonintphi[0][i]->GetNbinsX()+1),rebin[i],fittype);
 644     if(sigmafullmass[i]!=-2) fitter0.SetFixGaussianSigma(sigmafullmass[i]);
 645     ok=fitter0.MassFitter(kFALSE);
 646     if(ok){
 647       fitter0.DrawHere(cvnewptbins2->cd(i+1),nsigma,info);
 648       fitter0.Signal(nsigma,signal[0][i],errsgn[0][i]);
 649       fitter0.Background(nsigma,background[0][i],errbkg[0][i]);
 650       fitter0.Significance(nsigma,significance[0][i],errsignificance[0][i]);
 651
 652     }else{
 653       cvnewptbins2->cd(i+1);
 654       hmassnonintphi[0][i]->Draw();
 655     }
 656     cvnewptbins2->cd(i+1);
 657     pvbinpt->Draw();
 658
 659
 660     if(official) {
 661       cvnewptbins2->cd(i+1);
 662       TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
 663       txtsignif->SetBorderSize(0);
 664       txtsignif->SetFillStyle(0);
 665       txtsignif->SetTextColor(kGreen+3);
 666       txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,significance[0][i],errsignificance[0][i]));
 667
 668
 669       txtsignif->Draw();
 670       TH1F* htmp=(TH1F*)cvnewptbins2->cd(i+1)->FindObject("fhistoInvMass");
 671       htmp->SetYTitle(Form("Entries/%.3f GeV/c^{2}",htmp->GetBinWidth(3)));
 672       htmp->SetXTitle("D^{0} invariant mass (GeV/c^{2})");
 673     }
 674
 675     if(official) InOutPic(cvnewptbins2->cd(i+1),i);
 676
 677     //out-of-plane
 678     AliHFMassFitter fitter1(hmassnonintphi[1][i],hmassnonintphi[1][i]->GetBinLowEdge(1),hmassnonintphi[1][i]->GetBinLowEdge(hmassnonintphi[1][i]->GetNbinsX()+1),rebin[i],fittype);
 679     if(sigmafullmass[i]!=-2) fitter1.SetFixGaussianSigma(sigmafullmass[i]);
 680     ok=fitter1.MassFitter(kFALSE);
 681     if(ok){
 682       fitter1.DrawHere(cvnewptbins2->cd(nptbinsnew+i+1),nsigma,info);
 683       fitter1.Signal(nsigma,signal[1][i],errsgn[1][i]);
 684       fitter1.Background(nsigma,background[1][i],errbkg[1][i]);
 685       fitter1.Significance(nsigma,significance[1][i],errsignificance[1][i]);
 686     }else{
 687       cvnewptbins2->cd(nptbinsnew+i+1);
 688       hmassnonintphi[1][i]->Draw();
 689     }
 690     cvnewptbins2->cd(nptbinsnew+i+1);
 691     pvbinpt->Draw();
 692
 693     if(official) {
 694       cvnewptbins2->cd(nptbinsnew+i+1);
 695       TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
 696       txtsignif->SetBorderSize(0);
 697       txtsignif->SetFillStyle(0);
 698       txtsignif->SetTextColor(kGreen+3);
 699       txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,significance[1][i],errsignificance[1][i]));
 700
 701
 702       txtsignif->Draw();
 703       TH1F* htmp=(TH1F*)cvnewptbins2->cd(nptbinsnew+i+1)->FindObject("fhistoInvMass");
 704       htmp->SetYTitle(Form("Entries/%.3f GeV/c^{2}",htmp->GetBinWidth(3)));
 705       htmp->SetXTitle("D^{0} invariant mass (GeV/c^{2})");
 706     }
 707     if(official) InOutPic(cvnewptbins2->cd(nptbinsnew+i+1),i);
 708
 709   }
 710
 711   //other drawing details
 712   if(official){
 713     cvnewptbins2->cd(2);
 714     txtoff->Draw();
 715     gROOT->LoadMacro("/home/cbianchi/macros/ALICEPerformance.C");
 716     gROOT->ProcessLine(Form("ALICEPerformance((TVirtualPad*)%p,\"today\",\"br\")",cvnewptbins2->cd(1)));
 717   }
 718
 719
 720   //png mass histograms
 721   WriteCanvas(cvnewptbins2,Form("inout%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
 722   WriteCanvas(cvmasstotnewptbins,Form("totmass%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
 723   //graphs of yields
 724   TGraphErrors** gr2=new TGraphErrors*[2];
 725   TGraphErrors** grb2=new TGraphErrors*[2];
 726   TGraphErrors** grsgf2=new TGraphErrors*[2];
 727   TGraph** grerr2=new TGraph*[2];
 728   TGraph** grberr2=new TGraph*[2];
 729
 730   //canvas for these graphs
 731   TCanvas* cvsig2=new TCanvas("cvsig2","Signal vs pt (fixed sigma)");
 732   TCanvas* cvbkg2=new TCanvas("cvbkg2","Background vs pt (fixed sigma)");
 733   TCanvas* cvsgnf2=new TCanvas("cvsgnf2","Significance vs pt (fixed sigma)");
 734   TCanvas* cvsiger2=new TCanvas("cvsiger2","Error on Signal vs pt (fixed sigma)");
 735   TCanvas* cvbkger2=new TCanvas("cvbkger2","Error on Background vs pt (fixed sigma)");
 736
 737   for(Int_t i=0;i<2;i++) { //in/out
 738
 739     gr2[i]=new TGraphErrors(0);
 740     grb2[i]=new TGraphErrors(0);
 741     grsgf2[i]=new TGraphErrors(0);
 742     grerr2[i]=new TGraph(0);
 743     grberr2[i]=new TGraph(0);
 744
 745     for(Int_t j=0;j<nptbinsnew;j++){ //bins of pt
 746
 747       if(significance[i][j]>3.){
 748         gr2[i]->SetPoint(j,pt[j],signal[i][j]);
 749         gr2[i]->SetPointError(j,deltapt[j],errsgn[i][j]);
 750         grsgf2[i]->SetPoint(j,pt[j],significance[i][j]);
 751         grsgf2[i]->SetPointError(j,deltapt[j],errsignificance[i][j]);
 752         grerr2[i]->SetPoint(j,pt[j],errsgn[i][j]/signal[i][j]);
 753       }
 754       cout<<j<<" pt "<<pt[j]<<" pm "<<deltapt[j]<<", sgn "<<signal[i][j]<<" pm "<<errsgn[i][j]<<endl;
 755       grb2[i]->SetPoint(j,pt[j],background[i][j]);
 756       grb2[i]->SetPointError(j,deltapt[j],errbkg[i][j]);
 757       grberr2[i]->SetPoint(j,pt[j],errbkg[i][j]/background[i][j]);
 758     }
 759
 760     gr2[i]->SetTitle("Signal vs p_{t} (#sigma fixed)");
 761     gr2[i]->SetName(Form("sig%sfs", i==0 ? "inplane" : "outofplane"));
 762     gr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 763     gr2[i]->GetYaxis()->SetTitle("Signal");
 764     gr2[i]->SetMinimum(0);
 765     gr2[i]->SetMarkerStyle(20);
 766     gr2[i]->SetMarkerColor(i+1);
 767
 768     grb2[i]->SetTitle("Background vs p_{t} (#sigma fixed)");
 769     grb2[i]->SetName(Form("bkg%sfs", i==0 ? "inplane" : "outofplane"));
 770     grb2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 771     grb2[i]->GetYaxis()->SetTitle("Background");
 772     grb2[i]->SetMinimum(0);
 773     grb2[i]->SetMarkerStyle(20);
 774     grb2[i]->SetMarkerColor(i+1);
 775
 776     grsgf2[i]->SetTitle("Significance vs p_{t} (#sigma fixed);p_{t} (GeV/c);Significance");
 777     grsgf2[i]->SetName(Form("sgf%sfs", i==0 ? "inplane" : "outofplane"));
 778     grsgf2[i]->SetMinimum(0);
 779     grsgf2[i]->SetMarkerStyle(20);
 780     grsgf2[i]->SetMarkerColor(i+1);
 781
 782     grerr2[i]->SetTitle("#epsilon_{r} Signal (#sigma fixed)");
 783     grerr2[i]->SetName(Form("sgnerr%sfs", i==0 ? "inplane" : "outofplane"));
 784     grerr2[i]->GetYaxis()->SetTitle("#epsilon_{r} Signal");
 785     grerr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 786     grerr2[i]->SetMinimum(0);
 787     grerr2[i]->SetMarkerStyle(20);
 788     grerr2[i]->SetMarkerColor(i+1);
 789
 790     grberr2[i]->SetTitle("#epsilon_{r} Background (#sigma fixed)");
 791     grberr2[i]->SetName(Form("bkgerr%sfs", i==0 ? "inplane" : "outofplane"));
 792     grberr2[i]->GetYaxis()->SetTitle("#epsilon_{r} Background");
 793     grberr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 794     grberr2[i]->SetMinimum(0);
 795     grberr2[i]->SetMarkerStyle(20);
 796     grberr2[i]->SetMarkerColor(i+1);
 797
 798
 799     if(i==0) {
 800       cvsig2->cd();
 801       gr2[i]->Draw("AP");
 802
 803       cvbkg2->cd();
 804       grb2[i]->Draw("AP");
 805
 806       cvsgnf2->cd();
 807       grsgf2[i]->Draw("AP");
 808
 809       cvsiger2->cd();
 810       grerr2[i]->Draw("AP");
 811
 812       cvbkger2->cd();
 813       grberr2[i]->Draw("AP");
 814     }
 815     else {
 816       cvsig2->cd();
 817       gr2[i]->Draw("P");
 818
 819       cvbkg2->cd();
 820       grb2[i]->Draw("P");
 821
 822       cvsgnf2->cd();
 823       grsgf2[i]->Draw("P");
 824
 825       cvsiger2->cd();
 826       grerr2[i]->Draw("P");
 827
 828       cvbkger2->cd();
 829       grberr2[i]->Draw("P");
 830     }
 831
 832     //write in output file
 833     fout->cd();
 834     gr2[i]->Write();
 835     grb2[i]->Write();
 836     grerr2[i]->Write();
 837     grsgf2[i]->Write();
 838   }
 839
 840   cvsig2->cd();
 841   legsig->Draw();
 842   cvbkg2->cd();
 843   legsig->Draw();
 844   cvsgnf2->cd();
 845   legsig->Draw();
 846   cvsiger2->cd();
 847   legsig->Draw();
 848   cvbkger2->cd();
 849   legsig->Draw();
 850
 851   //write png
 852   WriteCanvas(cvsig2,textleg);
 853   WriteCanvas(cvbkg2,textleg);
 854   WriteCanvas(cvsgnf2,textleg);
 855   WriteCanvas(cvsiger2,textleg);
 856   WriteCanvas(cvbkger2,textleg);
 857
 858
 859   //invariant mass in bins of phi, pt integrated
 860
 861   Double_t sigvsphiptint[nphibins],errsigvsphiptint[nphibins],sgnfvsphiptint[nphibins],errsgnfvsphiptint[nphibins];
 862
 863   //canvas
 864   TCanvas* cvptint=new TCanvas("cvptint","Pt integrated mass plots",1600,600);
 865   cvptint->Divide(nphibins,1);
 866   TCanvas* cvphierr=new TCanvas("cvphierr","Error on Signal phi bins");
 867   TCanvas* cvchiphi=new TCanvas("cvchiphi","Chi as a function of phi");
 868   TCanvas* cvsgnvsphiptint=new TCanvas("cvsgnvsphiptint","Sigmal vs #Delta #phi");
 869
 870   //chi square and yields
 871   TH1F* hchivsphi=(TH1F*)hnphibins->Clone("hchivsphi");
 872   hchivsphi->SetTitle("#chi^{2} vs #phi;#phi (rad);#chi^{2}");
 873   hchivsphi->SetMarkerStyle(20);
 874   TGraphErrors* grphi=new TGraphErrors(0);
 875   TGraph* grphierr=new TGraph(0);
 876
 877   for (Int_t i=0;i<nphibins;i++){ //bins of phi
 878     cvptint->cd(i+1);
 879
 880     if(hmassptint[i]) {
 881       AliHFMassFitter fitter(hmassptint[i],hmassptint[i]->GetBinLowEdge(1),hmassptint[i]->GetBinLowEdge(hmassptint[i]->GetNbinsX()+1),rebin[0],fittype);
 882       Bool_t ok=fitter.MassFitter(kFALSE);
 883
 884       TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
 885       txtsignif->SetBorderSize(0);
 886       txtsignif->SetFillStyle(0);
 887       txtsignif->SetTextColor(kGreen+3);
 888
 889       if(ok) {
 890         fitter.DrawHere(cvptint->cd(i+1),nsigma,info);
 891         fitter.Signal(nsigma,sigvsphiptint[i],errsigvsphiptint[i]);
 892
 893         fitter.Significance(nsigma,sgnfvsphiptint[i],errsgnfvsphiptint[i]);
 894         txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,sgnfvsphiptint[i],errsgnfvsphiptint[i]));
 895
 896         //fill chi square
 897         hchivsphi->SetBinContent(i,fitter.GetChiSquare());
 898         //fill yields vs phi
 899         grphi->SetPoint(i,phibins[i],sigvsphiptint[i]);
 900         grphi->SetPointError(i,dx[i],errsigvsphiptint[i]);
 901         grphierr->SetPoint(i, phibins[i],errsigvsphiptint[i]/sigvsphiptint[i]);
 902       }
 903       else {
 904         //retry fit if it fails
 905         fitter.Reset();
 906         fitter.SetHisto(hmassptint[i]);
 907         fitter.SetRangeFit(hmassptint[i]->GetBinLowEdge(3),hmassptint[i]->GetBinLowEdge(hmassptint[i]->GetNbinsX()-2));
 908         fitter.RebinMass(rebin[0]*2);
 909         fitter.SetType(fittype,0);
 910         ok=fitter.MassFitter(kFALSE);
 911         if(ok){
 912           fitter.DrawHere(cvptint->cd(i+1),nsigma,info);
 913           TPaveText* pvnote=new TPaveText(0.1,0.7,0.8,0.9,"NDC");
 914           pvnote->SetBorderSize(0);
 915           pvnote->SetFillStyle(0);
 916           pvnote->AddText("Refitted w/ different range and rebin");
 917           cvptint->cd(i+1);
 918           pvnote->Draw();
 919
 920           fitter.Signal(nsigma,sigvsphiptint[i],errsigvsphiptint[i]);
 921
 922           //fill chi square
 923           hchivsphi->SetBinContent(i,fitter.GetChiSquare());
 924           //fill yields vs phi
 925           grphi->SetPoint(i,phibins[i],sigvsphiptint[i]);
 926           grphi->SetPointError(i,dx[i],errsigvsphiptint[i]);
 927           grphierr->SetPoint(i, phibins[i],errsigvsphiptint[i]/sigvsphiptint[i]);
 928         } else hmassptint[i]->Draw(); //fit didn't work at all
 929       }
 930       cvptint->cd(i+1);
 931
 932       if(official){
 933         txtsignif->Draw();
 934         //PhiBinPic(cvptint->cd(i+1),i,"tr");
 935       }
 936
 937     }
 938
 939   }
 940
 941   cvchiphi->cd();
 942   hchivsphi->Draw("ptext");
 943
 944   grphi->SetTitle("Signal vs #phi");
 945   grphi->SetName("gsigvsphi");
 946   grphi->GetXaxis()->SetTitle("#phi (rad)");
 947   grphi->GetYaxis()->SetTitle("Signal");
 948   grphi->SetMinimum(0);
 949   grphi->SetMarkerStyle(20);
 950   cvsgnvsphiptint->cd();
 951   grphi->Draw("AP");
 952
 953   grphierr->SetTitle("#epsilon_{r} Signal vs #phi");
 954   grphierr->SetName("gsigerrvsphi");
 955   grphierr->GetXaxis()->SetTitle("#phi (rad)");
 956   grphierr->GetYaxis()->SetTitle("#epsilon_{r} Signal");
 957   grphierr->SetMinimum(0);
 958   grphierr->SetMarkerStyle(20);
 959   cvphierr->cd();
 960   grphierr->Draw("AP");
 961
 962   //fit dN/dDeltaphi
 963   TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
 964   grphi->Fit(flowFunc);
 965
 966   TPaveText* pvv2=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
 967   pvv2->SetBorderSize(0);
 968   pvv2->SetFillStyle(0);
 969   pvv2->AddText(Form("v'_{2} = %.3f#pm%.3f",flowFunc->GetParameter(1),flowFunc->GetParError(1)));
 970
 971   cvsgnvsphiptint->cd();
 972   pvv2->Draw();
 973
 974   //write chi square, yields and errors in output file
 975   fout->cd();
 976   cvchiphi->Write();
 977   grphierr->Write();
 978   grphi->Write();
 979
 980   //other drawing details
 981   if(official){
 982     cvptint->cd(1);
 983     txtoff->Draw();
 984     gROOT->LoadMacro("/home/cbianchi/macros/ALICEPerformance.C");
 985     gROOT->ProcessLine(Form("ALICEPerformance((TVirtualPad*)%p,\"today\",\"br\")",cvptint->cd(1)));
 986   }
 987
 988   //write png mass histos yields and error
 989   WriteCanvas(cvptint,textleg);
 990   WriteCanvas(cvphierr,textleg);
 991   WriteCanvas(cvsgnvsphiptint,textleg);
 992
 993   //statistics
 994   TCanvas* cst=new TCanvas("cst","Stat");
 995   cst->SetGridy();
 996   cst->cd();
 997   hstat->Draw("htext0");
 998   WriteCanvas(cst,textleg);
 999
1000
1001   //save event plane info in the output file
1002   //fix to the latest version of the task!!!!!
1003   TH1F* hresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));
1004   /*
1005   TH2F* hMphis=(TH2F*)list->FindObject(Form("hMphicentr%d_%d",mincentr,mincentr+5));
1006   TH2F* hMc2phis=(TH2F*)list->FindObject(Form("hMc2phicentr%d_%d",mincentr,mincentr+5));
1007   */
1008   TString buildname=Form("centr%d_%d",mincentr,maxcentr);
1009   hresos->SetName(Form("hEvPlaneReso%s",buildname.Data()));
1010   /*
1011   hMphis->SetName(Form("hMphi%s",buildname.Data()));
1012   hMc2phis->SetName(Form("hMc2phi%s",buildname.Data()));
1013   */
1014   for(Int_t icentr=mincentr+5;icentr<=maxcentr;icentr=icentr+5){
1015     TH1F* hreso=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icentr,icentr+5));
1016     hresos->Add(hreso);
1017     // for(Int_t i=0;i<nptbins;i++){
1018     // TH2F* hMphi=(TH2F*)list->FindObject(Form("hMphicentr_pt%d%d_%d",icentr,icentr+5));
1019     // TH2F* hMc2phi=(TH2F*)list->FindObject(Form("hMc2phi_pt%dcentr%d_%d",icentr,icentr+5));
1020     // hMphis->Add(hMphi);
1021     // hMc2phis->Add(hMc2phi);
1022     // }
1023   }
1024
1025   fout->cd();
1026   hresos->Write();
1027   /*
1028   hMphis->Write();
1029   hMc2phis->Write();
1030   */
1031 }
1032
1033 void Dmesonsv2InOutAnisotropy(Int_t mincentr,Int_t maxcentr,Bool_t fixedsigma=kTRUE){
1034
1035   //Compute v2 with the method of the anisotropy between in-plane and out-of plane. Needs the output file from the previous method
1036   //Author: Francesco Prino prino@to.infn.it
1037
1038   gStyle->SetCanvasColor(0);
1039   gStyle->SetTitleFillColor(0);
1040   gStyle->SetStatColor(0);
1041   gStyle->SetOptTitle(0);
1042
1043   TString name=Form("centr%d-%d",mincentr,maxcentr),
1044     nameh=Form("centr%d_%d",mincentr,maxcentr);
1045
1046   TFile* fil=new TFile(Form("HistoInputv2Calc%s.root",name.Data()));
1047   if(!fil->IsOpen()){
1048     cout<<"Input file not found"<<endl;
1049     return;
1050   }
1051
1052
1053   //histogram resolution with subevents
1054   TH1F* hResolSubAB=(TH1F*)fil->Get(Form("hEvPlaneReso%s",nameh.Data()));
1055   hResolSubAB->GetXaxis()->SetTitle("cos[2(#psi_{A}-#psi_{B})]");
1056
1057   //graphs yield
1058   TGraphErrors* gsigin=(TGraphErrors*)fil->Get(Form("siginplane%s",fixedsigma ? "fs" : ""));
1059   TGraphErrors* gsigout=(TGraphErrors*)fil->Get(Form("sigoutofplane%s",fixedsigma ? "fs" : ""));
1060
1061   TCanvas* c1=new TCanvas("c1","EvPlaneResol");
1062   hResolSubAB->Draw();
1063   Double_t resolSub=TMath::Sqrt(hResolSubAB->GetMean());
1064   printf("Resolution on sub event  = %.4f\n",resolSub);
1065   Double_t chisub=FindChi(resolSub,1);
1066   printf("Chi Subevent             = %.4f\n",chisub);
1067   Double_t chifull=chisub*TMath::Sqrt(2);
1068   printf("Chi Full Event           = %.4f\n",chifull);
1069   Double_t resolFull=ComputeResol(resolSub,1);
1070   printf("Resolution on full event = %.4f\n",resolFull);
1071   TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
1072   tres->SetNDC();
1073   tres->Draw();
1074   c1->SaveAs(Form("EvPlaneResol%s.png",nameh.Data()));
1075
1076   TCanvas* c2=new TCanvas("c2",Form("Yields %s",fixedsigma ? "fixed sigma" : ""));
1077   c2->cd();
1078   gsigin->Draw("AP");
1079   gsigout->Draw("Psame");
1080
1081   Int_t nPtBins=gsigin->GetN();
1082
1083   //graphs results
1084   TGraphErrors *gAnis=new TGraphErrors(0);
1085   gAnis->SetName(Form("gAnisotropy%s",nameh.Data()));
1086   gAnis->SetTitle("(Nin-Nout)/(Nin+Nout)");
1087
1088   TGraphErrors *gv2obs=new TGraphErrors(0);
1089   gv2obs->SetName(Form("gv2obs%s",nameh.Data()));
1090   gv2obs->SetTitle("v2obs");
1091
1092   TGraphErrors *gv2=new TGraphErrors(0);
1093   gv2->SetName(Form("gv2%s",nameh.Data()));
1094   gv2->SetTitle("v2");
1095
1096   TGraph *gv2err=new TGraph(0);
1097   gv2err->SetName(Form("gv2err%s",nameh.Data()));
1098
1099   for(Int_t iBin=0; iBin<nPtBins; iBin++){
1100     Double_t pt,nIn,ept,enIn;
1101     Double_t nOut,enOut;
1102     gsigin->GetPoint(iBin,pt,nIn);
1103     ept=gsigin->GetErrorX(iBin);
1104     enIn=gsigin->GetErrorY(iBin);
1105     gsigout->GetPoint(iBin,pt,nOut);
1106     enOut=gsigout->GetErrorY(iBin);
1107     if(nIn<1e-6 || nOut<1e-6) continue;
1108
1109     Double_t anis=(nIn-nOut)/(nIn+nOut);
1110     Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
1111     gAnis->SetPoint(iBin,pt,anis);
1112     gAnis->SetPointError(iBin,ept,eAnis);
1113
1114     Double_t v2obs=anis*TMath::Pi()/4.;
1115     Double_t ev2obs=eAnis*TMath::Pi()/4.;
1116     gv2obs->SetPoint(iBin,pt,v2obs);
1117     gv2obs->SetPointError(iBin,ept,ev2obs);
1118
1119     Double_t v2=v2obs/resolFull;
1120     Double_t ev2=ev2obs/resolFull;
1121     gv2->SetPoint(iBin,pt,v2);
1122     gv2->SetPointError(iBin,ept,ev2);
1123
1124     gv2err->SetPoint(iBin,pt,ev2/TMath::Abs(v2));
1125   }
1126
1127   TCanvas* c3=new TCanvas("c3","Anisotropy");
1128   gAnis->SetMarkerStyle(20);
1129   gAnis->GetXaxis()->SetTitle("p_{t} (GeV/c)");
1130   gAnis->GetYaxis()->SetTitle("(N_{IN-PLANE}-N_{OUT-OF-PLANE})/(N_{IN-PLANE}+N_{OUT-OF-PLANE})");
1131   gAnis->GetYaxis()->SetTitleOffset(1.1);
1132   gAnis->Draw("AP");
1133   gAnis->Draw("Psame");
1134   c3->SaveAs(Form("Anisotropy%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
1135
1136   TCanvas* c4=new TCanvas("c4","v2obs");
1137   gv2obs->SetMarkerStyle(20);
1138   gv2obs->GetXaxis()->SetTitle("p_{t} (GeV/c)");
1139   gv2obs->GetYaxis()->SetTitle("v_{2}^{obs}");
1140   gv2obs->GetYaxis()->SetTitleOffset(1.1);
1141   gv2obs->Draw("AP");
1142   gv2obs->Draw("Psame");
1143   c4->SaveAs(Form("v2obsD%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
1144
1145   TCanvas* c5=new TCanvas("c5","v2");
1146   gv2->SetMarkerStyle(20);
1147   gv2->GetXaxis()->SetTitle("p_{t} (GeV/c)");
1148   gv2->GetYaxis()->SetTitle("v_{2}");
1149   gv2->GetYaxis()->SetTitleOffset(1.1);
1150   gv2->Draw("AP");
1151   gv2->Draw("Psame");
1152   c5->SaveAs(Form("v2D%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
1153
1154   TCanvas* c6=new TCanvas("c6","error v2");
1155   gv2err->SetTitle("Relative error v2;p_{t} (GeV/c);#epsilon_{r} v2");
1156   gv2err->SetMarkerStyle(20);
1157   gv2err->GetYaxis()->SetTitleOffset(1.1);
1158   c6->cd();
1159   gv2err->Draw("AP");
1160   c6->SaveAs(Form("Errv2%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
1161
1162   TFile* fout=new TFile(Form("Outputv2%s.root",name.Data()),"recreate");
1163
1164   fout->cd();
1165   gAnis->Write();
1166   gv2obs->Write();
1167   gv2->Write();
1168   gv2err->Write();
1169   c1->Write(); //canvas resolution
1170 }
1171
1172 //methods needed inside Dmesonsv2InOutAnisotropy
1173 //Author: Francesco Prino prino@to.infn.it
1174
1175 Double_t ComputeResol(Double_t resSub, Int_t k){
1176   Double_t chisub=FindChi(resSub,k);
1177   Double_t chifull=chisub*TMath::Sqrt(2);
1178   if(k==1) return ResolK1(chifull);
1179   else if(k==2) return Pol(chifull,2);
1180   else{
1181     printf("k should be <=2\n");
1182     return 1.;
1183   }
1184 }
1185
1186 Double_t FindChi(Double_t res, Int_t k){
1187   Double_t x1=0;
1188   Double_t x2=15;
1189   Double_t y1,y2;
1190   if(k==1){
1191     y1=ResolK1(x1)-res;
1192     y2=ResolK1(x2)-res;
1193   }
1194   else if(k==2){
1195     y1=Pol(x1,2)-res;
1196     y2=Pol(x2,2)-res;
1197   }
1198   else return -1;
1199
1200   if(y1*y2>0.) return -1;
1201   if(y1==0.) return y1;
1202   if(y2==0.) return y2;
1203   Double_t xmed,ymed;
1204   Int_t jiter=0;
1205   while((x2-x1)>0.0001){
1206     xmed=0.5*(x1+x2);
1207     if(k==1){
1208       y1=ResolK1(x1)-res;
1209       y2=ResolK1(x2)-res;
1210       ymed=ResolK1(xmed)-res;
1211     }
1212     else if(k==2){
1213       y1=Pol(x1,2)-res;
1214       y2=Pol(x2,2)-res;
1215       ymed=Pol(xmed,2)-res;
1216     }
1217     else return -1;
1218     if((y1*ymed)<0) x2=xmed;
1219     if((y2*ymed)<0)x1=xmed;
1220     if(ymed==0) return xmed;
1221     jiter++;
1222   }
1223   return 0.5*(x1+x2);
1224 }
1225
1226 Double_t Pol(Double_t x, Int_t k){
1227   Double_t c[5];
1228   if(k==1){
1229     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
1230   }
1231   else if(k==2){
1232     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
1233   }
1234   else return -1;
1235   return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
1236 }
1237
1238 Double_t ResolK1(Double_t x){
1239   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
1240 }
1241
1242
1243 void DrawEventPlane(Int_t mincentr=0,Int_t maxcentr=0,TString partname="D0",TString path="./"){
1244
1245  //draw the histograms correlated to the event plane
1246
1247   gStyle->SetCanvasColor(0);
1248   gStyle->SetTitleFillColor(0);
1249   gStyle->SetStatColor(0);
1250   //gStyle->SetPalette(1);
1251   gStyle->SetOptFit(1);
1252
1253   TString listname="coutputv2";
1254   TList* list;
1255   TH1F * hstat;
1256   AliRDHFCuts* cutobj;
1257
1258   Bool_t isRead=ReadFile(list,hstat,cutobj,listname,partname,path);
1259   if(!isRead) return;
1260   if(!list || !hstat || !cutobj){
1261     cout<<":-( null pointers..."<<endl;
1262     return;
1263   }
1264
1265   if(!(mincentr==0 && maxcentr==0)){ //draw the total in mincentr-maxcentr
1266     TString suffixcentr=Form("centr%d_%d",mincentr,maxcentr);
1267     TH1F* hevpls=(TH1F*)list->FindObject(Form("hEvPlanecentr%d_%d",mincentr,mincentr+5));
1268     hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
1269     hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
1270     TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));
1271     hevplresos->SetName(Form("hEvPlaneReso%s",suffixcentr.Data()));
1272     hevplresos->SetTitle(Form("Event Plane Resolution %s",suffixcentr.Data()));
1273
1274     for(Int_t icentr=mincentr+5;icentr<=maxcentr;icentr=icentr+5){
1275       TH1F* h=(TH1F*)list->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+5));
1276       if(h)hevpls->Add(h);
1277       else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
1278       h=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icentr,icentr+5));
1279       if(h)hevplresos->Add(h);
1280       else cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+5<<endl;
1281     }
1282
1283     TCanvas* cvtotevpl=new TCanvas("cvtotevpl",Form("Ev plane %s",suffixcentr.Data()));
1284     cvtotevpl->cd();
1285     hevpls->Draw();
1286     hevpls->Fit("pol0");
1287
1288     TCanvas* cvtotevplreso=new TCanvas("cvtotevplreso",Form("Ev plane Resolution %s",suffixcentr.Data()));
1289     cvtotevplreso->cd();
1290     hevplresos->Draw();
1291     Double_t resolSub=TMath::Sqrt(hevplresos->GetMean());
1292     Double_t resolFull=ComputeResol(resolSub,1);
1293
1294     TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
1295     pvreso->SetBorderSize(0);
1296     pvreso->SetFillStyle(0);
1297     pvreso->AddText(Form("Resolution on full event = %.4f\n",resolFull));
1298     pvreso->Draw();
1299
1300     TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",mincentr,maxcentr),"recreate");
1301     fout->cd();
1302     hevpls->Write();
1303     hevplresos->Write();
1304   }
1305   else{ //draw all in 5% centrality bins
1306     for(Int_t i=0;i<list->GetEntries();i++){
1307       TH1F* h=(TH1F*)list->At(i);
1308       if(!h){
1309         cout<<"Histogram "<<i<<" not found"<<endl;
1310         continue;
1311       }
1312       TString hname=h->GetName();
1313       if(hname.Contains("EvPlane")){
1314         TCanvas* cv=new TCanvas(Form("cv%s",hname.Data()),hname.Data());
1315         cv->cd();
1316         h->Draw();
1317
1318         if(!hname.Contains("Reso")){
1319           h->Fit("pol0");
1320
1321         }else{
1322           Double_t resolSub=TMath::Sqrt(h->GetMean());
1323           Double_t resolFull=ComputeResol(resolSub,1);
1324
1325           TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
1326           pvreso->SetBorderSize(0);
1327           pvreso->SetFillStyle(0);
1328           pvreso->AddText(Form("Resolution on full event = %.4f\n",resolFull));
1329           pvreso->Draw();
1330
1331         }
1332       }
1333     }
1334   }
1335
1336 }
1337
1338
1339 void DmesonsFlowAnalysisvsCentrality(Int_t mincentr, Int_t maxcentr, Int_t step=10){
1340   //needs outputs from the other methods HistoInputv2Calccentr%d-%d.root,EvPlanecentr%d-%d.root, Outputv2centr%d-%d.root
1341
1342   //min(max)centr= min(max) value of centrality among those files - note that if there are "holes" they are skipped
1343   //step= step in centrality between mincentr and maxcentr - default 10%
1344
1345   gStyle->SetCanvasColor(0);
1346   gStyle->SetTitleFillColor(0);
1347   gStyle->SetStatColor(0);
1348   gStyle->SetOptStat(0);
1349   gStyle->SetOptTitle(0);
1350
1351   //number of files expected
1352   Int_t n=(maxcentr-mincentr)/step;
1353   cout<<"Trying to read "<<n<<" files"<<endl;
1354
1355   TString filename="",gname="";
1356   Int_t icentr=0;
1357
1358   //v2
1359   TCanvas* cvv2=new TCanvas("cvv2","v_{2}");
1360   TCanvas* cvv2obs=new TCanvas("cvv2obs","v_{2}'");
1361
1362   //Event plane distribution and resolution
1363   TCanvas* cvevplres=new TCanvas("cvevplres","Resolution");
1364   cvevplres->SetLogy();
1365   TH1F* hresolFull=new TH1F("hresolFull","Resolution vs centrality;centrality;resolution",n,mincentr,maxcentr);
1366   TCanvas* cvres=new TCanvas("cvres", "Resolution vs centrality");
1367   TCanvas* cvevplane=new TCanvas("cvevplane", "Event plane vs centrality");
1368
1369   //v2 in phi bins
1370   TCanvas* cvfitvsphi=new TCanvas("cvfitvsphi","Fit flow vs phi");
1371   TH1F* hv2fcorr=new TH1F("hv2fcorr","v2 from fit corrected;centrality;v_{2}",n,mincentr,maxcentr);
1372   hv2fcorr->SetMarkerStyle(30);
1373   hv2fcorr->SetMarkerSize(2);
1374   TH1F* hv2corr=new TH1F("hv2corr","v2 from asym corrected;centrality;v_{2}",n,mincentr,maxcentr);
1375   hv2corr->SetMarkerStyle(29);
1376   hv2corr->SetMarkerSize(1.7);
1377   hv2corr->SetMarkerColor(2);
1378   hv2corr->SetLineColor(2);
1379
1380   TCanvas* cvv2cor=new TCanvas("v2corr", "v2 from asymmetry corrected vs centrality");
1381
1382   TLegend* leg=new TLegend(0.7,0.6,0.9,0.8,"Centrality %");
1383   leg->SetBorderSize(0);
1384   leg->SetFillStyle(0);
1385   TLegend* leg2=new TLegend(0.5,0.6,0.9,0.9,"");
1386   leg2->SetBorderSize(0);
1387   leg2->SetFillStyle(0);
1388   TLegend* legresoval=new TLegend(0.2,0.6,0.5,0.9,"Full resolution");
1389   legresoval->SetBorderSize(0);
1390   legresoval->SetFillStyle(0);
1391
1392   Int_t nread=0;
1393   for(Int_t i=0;i<n;i++){
1394     icentr=mincentr+i*step;
1395     filename=Form("Outputv2centr%d-%d.root",icentr,icentr+step);
1396     TFile* f=new TFile(filename);
1397     if(!f->IsOpen()){
1398       cout<<filename.Data()<<" not found"<<endl;
1399       continue;
1400     }
1401     nread++;
1402     gname=Form("gv2centr%d_%d",icentr,icentr+step);
1403     TGraphErrors* grv2=(TGraphErrors*)f->Get(gname);
1404     if(!grv2){
1405       cout<<gname<<" not found in "<<filename<<endl;
1406       continue;
1407     }
1408     grv2->GetYaxis()->SetRangeUser(-0.2,0.4);
1409     grv2->SetMarkerColor(i+1);
1410     grv2->SetLineColor(i+1);
1411     cvv2->cd();
1412     if(nread==1)grv2->Draw("APL");
1413     else grv2->Draw("PL");
1414
1415     leg->AddEntry(grv2,Form("%d-%d",icentr,icentr+step),"p");
1416
1417     gname=Form("gv2obscentr%d_%d",icentr,icentr+step);
1418     TGraphErrors* grv2obs=(TGraphErrors*)f->Get(gname);
1419     if(!grv2obs){
1420       cout<<gname<<" not found in "<<filename<<endl;
1421       continue;
1422     }
1423     grv2obs->GetYaxis()->SetRangeUser(-0.2,0.4);
1424     grv2obs->SetMarkerColor(i+1);
1425     grv2obs->SetLineColor(i+1);
1426     cvv2obs->cd();
1427     if(nread==1)grv2obs->Draw("APL");
1428     else grv2obs->Draw("PL");
1429
1430     //draw ev plane reso
1431     filename=Form("EvPlanecentr%d-%d.root",icentr,icentr+step);
1432     TFile* f2=new TFile(filename);
1433     if(!f2->IsOpen()){
1434       cout<<filename.Data()<<" not found"<<endl;
1435       continue;
1436     }
1437     TString hname=Form("hEvPlaneResocentr%d_%d",icentr,icentr+step);
1438     TH1F* hevplres=(TH1F*)f2->Get(hname);
1439     if(!hevplres){
1440       cout<<hname<<" not found in "<<filename<<endl;
1441       continue;
1442     }
1443     cvevplres->cd();
1444     hevplres->SetLineColor(i+1);
1445     if(nread==1) hevplres->Draw();
1446     else hevplres->Draw("sames");
1447
1448     Double_t resolSub=TMath::Sqrt(hevplres->GetMean());
1449     Double_t resolFull=ComputeResol(resolSub,1);
1450     hresolFull->SetBinContent(i+1,resolFull);
1451     TLegendEntry* lentry=legresoval->AddEntry(hevplres,Form("%.2f",resolFull),"");
1452     lentry->SetTextColor(i+1);
1453     cvevplres->cd();
1454     legresoval->Draw();
1455
1456     hname=Form("hEvPlanecentr%d_%d",icentr,icentr+step);
1457     TH1F* hevpl=(TH1F*)f2->Get(hname);
1458     if(!hevpl){
1459       cout<<hname<<" not found in "<<filename<<endl;
1460       continue;
1461     }
1462     hevpl->Scale(1./hevpl->Integral());
1463     cvevplane->cd();
1464     hevpl->SetLineColor(i+1);
1465     if(nread==1) hevpl->Draw();
1466     else hevpl->Draw("sames");
1467
1468
1469     filename=Form("HistoInputv2Calccentr%d-%d.root",icentr,icentr+step);
1470     TFile* f3=new TFile(filename);
1471     if(!f3->IsOpen()){
1472       cout<<filename.Data()<<" not found"<<endl;
1473       continue;
1474     }
1475     gname="gsigvsphi";
1476     TGraphErrors* grdNdphi=(TGraphErrors*)f3->Get(gname);
1477     if(!grdNdphi){
1478       cout<<gname<<" not found in "<<filename<<endl;
1479       continue;
1480     }
1481     TF1* ffit=(TF1*)grdNdphi->FindObject("flow");
1482     ffit->SetName(Form("flow%d_%d",icentr,icentr+step));
1483     ffit->SetLineColor(i+1);
1484     // Double_t xmin,xmax;
1485     // ffit->GetRange(xmin,xmax);
1486     //ffit->Scale(1./ffit->Integral(xmin,xmax));
1487     cvfitvsphi->cd();
1488     if(nread==1){
1489       ffit->SetMinimum(0);
1490       ffit->Draw();
1491     }
1492     else ffit->Draw("sames");
1493
1494     Double_t v2obsf=ffit->GetParameter(1);
1495     Double_t v2obsferr=ffit->GetParError(1);
1496     Double_t v2fcorr=v2obsf/resolFull;
1497     Double_t v2fcorrerr=v2obsferr/resolFull;
1498     hv2fcorr->SetBinContent(i+1,v2fcorr);
1499     hv2fcorr->SetBinError(i+1,v2fcorrerr);
1500
1501     gname="siginplane";
1502     TGraphErrors* gsigin=(TGraphErrors*)f3->Get(gname);
1503     gname="sigoutofplane";
1504     TGraphErrors* gsigout=(TGraphErrors*)f3->Get(gname);
1505     if(!gsigin || !gsigout){
1506       cout<<"graphs in/out not found"<<endl;
1507       continue;
1508     }
1509     Int_t nbins=((TParameter<float>*)f3->Get("nptbins"))->GetVal();
1510     Double_t countin=0,countout=0,e2in=0,e2out=0;
1511     for(Int_t k=0;k<nbins;k++){
1512       Double_t x,y,ey;
1513       gsigin->GetPoint(k,x,y);
1514       ey=gsigin->GetErrorY(k);
1515       e2in+=ey*ey;
1516       countin+=y;
1517       gsigout->GetPoint(k,x,y);
1518       ey=gsigin->GetErrorY(k);
1519       e2out+=ey*ey;
1520       countout+=y;
1521     }
1522     if(countin<1e-6 || countout<1e-6) continue;
1523     Double_t ein=TMath::Sqrt(e2in),eout=TMath::Sqrt(e2out);
1524     Double_t anis=(countin-countout)/(countin+countout);
1525     Double_t eAnis=TMath::Sqrt(ein*ein+eout*eout)/(countin+countout);
1526     Double_t v2obs=anis*TMath::Pi()/4.;
1527     Double_t ev2obs=eAnis*TMath::Pi()/4.;
1528     Double_t v2corr=v2obs/resolFull;
1529     Double_t v2correrr=ev2obs/resolFull;
1530     hv2corr->SetBinContent(i+1,v2corr);
1531     hv2corr->SetBinError(i+1,v2correrr);
1532
1533   }
1534
1535   cvv2->cd();
1536   leg->Draw();
1537   cvv2obs->cd();
1538   leg->Draw();
1539   cvevplres->cd();
1540   leg->Draw();
1541   cvevplane->cd();
1542   leg->Draw();
1543
1544   cvres->cd();
1545   hresolFull->Draw();
1546
1547   cvfitvsphi->cd();
1548   leg->Draw();
1549
1550   cvv2cor->cd();
1551   hv2corr->SetLineWidth(2);
1552   hv2fcorr->SetLineWidth(2);
1553   hv2corr->GetYaxis()->SetRangeUser(-0.2,0.4);
1554   hv2corr->Draw();
1555   hv2fcorr->Draw("sames");
1556   leg2->AddEntry(hv2fcorr,"v_{2} from dN/d#phi","l");
1557   leg2->AddEntry(hv2corr,"v_{2} from asymmetry","l");
1558   leg2->Draw();
1559
1560   cvv2->SaveAs(Form("v2Dcentr%d-%d.png",mincentr, maxcentr));
1561   cvv2obs->SaveAs(Form("v2obsDcentr%d-%d.png",mincentr, maxcentr));
1562   cvevplres->SaveAs(Form("evplaneresocentr%d-%d.png",mincentr, maxcentr));
1563   cvevplane->SaveAs(Form("evplanedistrcentr%d-%d.png",mincentr, maxcentr));
1564   cvres->SaveAs(Form("evplresovscentr%d-%d.png",mincentr, maxcentr));
1565   cvfitvsphi->SaveAs(Form("flowvsphicentr%d-%d.png",mincentr, maxcentr));
1566   cvv2cor->SaveAs(Form("cmpv2centr%d-%d.png",mincentr, maxcentr));
1567
1568
1569 }
1570
1571 void ReflectdNdphiPoints(Int_t mincentr,Int_t maxcentr){
1572
1573   gStyle->SetCanvasColor(0);
1574   gStyle->SetTitleFillColor(0);
1575   gStyle->SetOptTitle(0);
1576   gStyle->SetStatColor(0);
1577   gStyle->SetOptStat(0);
1578
1579   TString textleg=Form("centr%d-%d",mincentr,maxcentr);
1580   TFile* fin=new TFile(Form("HistoInputv2Calc%s.root",textleg.Data()));
1581   if(!fin->IsOpen()){
1582     printf("HistoInputv2Calc%s.root not found",textleg.Data());
1583     return;
1584   }
1585   TGraphErrors* grphi=(TGraphErrors*)fin->Get("gsigvsphi");
1586   if(!grphi){
1587     cout<<"gsigvsphi not found in "<<fin->GetName()<<endl;
1588     return;
1589   }
1590
1591   TGraphErrors* grphidouble=new TGraphErrors(0);
1592   Int_t n=grphi->GetN();
1593   //Int_t doublen=n*2;
1594   Double_t xn=0,yn=0;
1595   grphi->GetPoint(n-1,xn,yn);
1596   xn+=grphi->GetErrorX(n-1);
1597   cout<<"xn = "<<xn<<endl;
1598   for(Int_t i=0;i<n;i++){
1599     Double_t x,y;
1600     grphi->GetPoint(i,x,y);
1601     //x+=xn;
1602     x=(-1)*x;
1603     grphidouble->SetPoint(i,x,y);
1604     grphidouble->SetPointError(i,grphi->GetErrorX(i),grphi->GetErrorY(i));
1605   }
1606   grphidouble->SetMarkerStyle(24);
1607   grphidouble->SetName(Form("%sReflected",grphi->GetName()));
1608   grphidouble->SetTitle(grphi->GetTitle());
1609   grphidouble->GetXaxis()->SetRangeUser(0,6.5);
1610
1611   TMultiGraph* grtot=new TMultiGraph();
1612   grtot->Add(grphi,"P");
1613   grtot->Add(grphidouble,"P");
1614   grtot->SetNameTitle("grphi2pi","dN/d#Delta#phi;#Delta#phi (rad);dN/d#Delta#phi");
1615   TPaveText* pv=new TPaveText(0.3,0.7,0.7,0.9,"NDC");
1616   pv->SetBorderSize(0);
1617   pv->SetFillStyle(0);
1618
1619   TPaveText* pvcentr=new TPaveText(0.1,0.7,0.5,0.8,"NDC");
1620   pvcentr->SetBorderSize(0);
1621   pvcentr->SetFillStyle(0);
1622   pvcentr->AddText(Form("Centrality %d-%d",mincentr,maxcentr));
1623
1624   TPaveText* pvwp=new TPaveText(0.5,0.1,0.9,0.3,"NDC");
1625   pvwp->SetBorderSize(0);
1626   pvwp->SetFillStyle(0);
1627   pvwp->SetTextColor(kRed);
1628   pvwp->SetTextFont(65);
1629   pvwp->AddText("ALICE Work in progress");
1630
1631   pv->AddText("Open points are reflected");
1632   TCanvas *cvphidouble=new TCanvas("cvphidouble","dN/dphi in 2*pi");
1633   cvphidouble->cd();
1634   grtot->Draw("AP");
1635   pv->Draw();
1636
1637   TCanvas *cvphi=new TCanvas("cvphi","dN/dphi");
1638   cvphi->cd();
1639   grphi->SetNameTitle(grphi->GetName(),"dN/d#Delta#phi;#Delta#phi (rad);dN/d#Delta#phi");
1640   grphi->Draw("AP");
1641   pvcentr->Draw();
1642   pvwp->Draw();
1643 }